Simulating LTI Systems in Matlab
This tutorial will assume you have created state-space matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ which characterize the system modeled by
$$ \mathbf{\dot{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t) $$ $$ \mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}(t) $$and have stored these matrices in MATLAB variables A
, B
, C
, and D
, respectively. The time-domain response of the system can then be solved numerically in MATLAB using two functions, ss()
and lsim()
. The first few lines of the help documentation for each is given below
>> help ss
SYS = ss(A,B,C,D) creates an object SYS representing the continuous-time state-space model
dx/dt = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)
>>help lsim
lsim(SYS,U,T) plots the time response of the dynamic system SYS to the input signal described by U and T.
The time vector T is expressed in the time units of SYS and consists of regularly spaced time samples.
The matrix U has as many columns as inputs in SYS and its i-th row specifies the input value at time T(i). For example,
t = 0:0.01:5; u = sin(t); lsim(sys,u,t)
simulates the response of a single-input model SYS to the input u(t)=sin(t) during 5 time units.
For a fully LTI system, where all state-space matrices are constant, the solution at time $t>0$ can be solved numerically by simply integrating $\mathbf{\dot{x}}(t)$ over time, beginning from some initial condition. For example, for the problem of homework #2, one way to setup the simulation is as follows
sys = ss(A,B,C,D); %combines matrices to generate a system
nPeriods = 200; %number of periods to simulate
pts_per_period = 10000; %number of time-points per period
npts = nPeriods*pts_per_period; %total number of time-points in simulation
t = linspace(0, nPeriods*Ts, npts); %time vector for simulation
vp_Ts = [-Vg*ones(1,pts_per_period/2),...
Vg*ones(1,pts_per_period/2)]; %creates one period of squarewave vp(t)
vp = repmat(vp_Ts, 1, nPeriods); %repeats the single period of vp(t) nPeriods times
vs = vp/Vg*V; %creates vs(t). This just re-uses vp(t) and changes the amplitude
vp = circshift(vp, round(tphi/Ts*pts_per_period)); %adds a phase-shift of tphi to vp
u = cat(1,vp,vs); %concatenate vp and vs together for the u(t) vector needed in lsim
[Y,T,X] = lsim(sys,u,t,zeros(nstates,1)); %simulates the system with the input u(t), beginning at zero initial conditions