lsim

假设有以下传递函数。

sys = tf(3,[1 2 3])

sys =
 
        3
  -------------
  s^2 + 2 s + 3
 
Continuous-time transfer function.
Model Properties

要计算此系统对任意输入信号的响应，需向 lsim 提供一个包含您要计算响应的时间的时间向量 t，以及一个包含相应信号值的向量 u。例如，绘制此系统对斜坡阶跃信号的响应，该信号在时间 t = 0 时从 0 开始，在 t = 1 时从 0 斜坡上升到 t = 2 时的 1，然后保持稳定在 1。定义 t 并计算 u 的值。

t = 0:0.04:8;  % 201 points
u = max(0,min(t-1,1));

不使用任何输出参量调用 lsim 来绘制此系统对该信号的响应。

lsim(sys,u,t)
grid on

该图以灰色显示应用的输入 (u,t)，以蓝色显示系统响应。

使用输出参量调用 lsim 来获取仿真响应数据。

y = lsim(sys,u,t);
size(y)

ans = 1×2

   201     1

向量 y 包含在 t 中相应时间点的仿真响应。

使用 gensig (Control System Toolbox) 创建周期性输入信号（如正弦波和方波）以用于 lsim。仿真以下 SISO 状态空间模型对方波的响应。

A = [-3 -1.5; 5 0];
B = [1; 0];
C = [0.5 1.5];
D = 0;
sys = ss(A,B,C,D);

对于此示例，创建一个周期为 10 秒、持续时间为 20 秒的方波。

[u,t] = gensig("square",10,20);

gensig 返回时间步的向量 t 和包含输入信号相应值的向量 u。（如果您没有为 t 指定采样时间，则 gensig 会每个周期生成 64 个采样。）对 lsim 使用这两个向量，并绘制系统响应。

lsim(sys,u,t)
grid on

该图以灰色显示施加的方波，以蓝色显示系统响应。使用输出参量调用 lsim 来获取 t 中每个时间点的响应值。

[y,~] = lsim(sys,u,t);

对离散时间系统的响应进行仿真时，时间向量 t 必须采用 Ti:dT:Tf 的形式，其中 dT 是模型的采样时间。仿真以下离散时间传递函数对斜坡阶跃输入的响应。

sys = tf([0.06 0.05],[1 -1.56 0.67],0.05);

此传递函数的采样时间为 0.05 秒。使用相同的采样时间生成时间向量 t 和斜坡阶跃信号 u。

t = 0:0.05:4;  
u = max(0,min(t-1,1));

绘制系统响应。

lsim(sys,u,t)

要仿真离散时间系统对周期性输入信号的响应，请将相同的采样时间用于 gensig 来生成输入。例如，仿真系统对周期为 1 秒、持续时间为 4 秒的正弦波的响应。

[u,t] = gensig("sine",1,4,0.05);

绘制系统响应。

lsim(sys,u,t)

使用 lsim，您可以在同一轴上绘制多个动态系统的仿真响应。例如，将系统的闭环响应与 PI 控制器和 PID 控制器的闭环响应进行比较。创建系统的传递函数并调节控制器。

H = tf(4,[1 10 25]);
C1 = pidtune(H,'PI');
C2 = pidtune(H,'PID');

形成闭环系统。

sys1 = feedback(H*C1,1);
sys2 = feedback(H*C2,1);

绘制两个系统对周期为 4 秒的方波的响应。

[u,t] = gensig("square",4,12);
lsim(sys1,sys2,u,t)
grid on
legend("PI","PID")

默认情况下，lsim 会为您绘制的每个系统选择不同的颜色。您可以使用 LineSpec 输入参量指定颜色和线型。

 lsim(sys1,"r--",sys2,"b",u,t)
 grid on
 legend("PI","PID")

第一个 LineSpec "r--" 为 PI 控制器的响应指定红色虚线。第二个 LineSpec "b" 为 PID 控制器的响应指定蓝色实线。图例反映了所指定的颜色和线型。如需更多绘图自定义选项，请使用 lsimplot。

在 MIMO 系统中，在每个时间步 t 处，输入 u(t) 是长度等于输入数的向量。要使用 lsim，您需要将 u 指定为维度为 Nt×Nu 的矩阵，其中 Nu 是系统输入的数目，Nt 是 t 的长度。换句话说，u 的每一列都是施加于相应系统输入端的输入信号。例如，要在 201 个时间步对一个具有四个输入的系统进行仿真，需要提供一个由四列和 201 行组成的矩阵作为 u，其中每一行 u(i,:) 是第 i 个时间步时的输入值向量；每一列 u(:,j) 是在第 j 个输入端处施加的信号。

类似地，由 lsim 计算得到的输出 y(t) 是一个矩阵，其列表示每个系统输出端处的信号。当您使用 lsim 绘制仿真响应时，lsim 会为每个输出提供单独的坐标区，表示每个输出通道中的对在所有输入端施加的输入 u(t) 的系统响应。

假设有具有以下状态空间矩阵的双输入、三输出状态空间模型。

A = [-1.5  -0.2   1.0;
     -0.2  -1.7   0.6;
      1.0   0.6  -1.4];
  
B = [ 1.5  0.6;
     -1.8  1.0;
      0    0  ];

C = [ 0    -0.5 -0.1;
      0.35 -0.1 -0.15
      0.65  0    0.6];
  
D = [ 0.5  0;
      0.05 0.75
      0    0];

sys = ss(A,B,C,D);

绘制 sys 对周期为 4 秒的方波（施加到第一个输入 sys）和每 3 秒施加一次的冲激（施加到第二个输入）的响应。为此，使用 gensig 创建表示该方波和该冲激信号的列向量。然后将这些列堆叠组合成一个输入矩阵。为确保两个信号具有相同的采样数，需指定相同的结束时间和采样时间。

Tf = 10;
Ts = 0.1;
[uSq,t] = gensig("square",4,Tf,Ts);
[uP,~] = gensig("pulse",3,Tf,Ts);
u = [uSq uP];
lsim(sys,u,t)

每个轴显示三个系统输出端之一对施加在所有输入端的信号 u 的响应。每个图还以灰色显示所有输入信号。

默认情况下，lsim 在仿真模型时假设所有状态在仿真开始时均为零。在对状态空间模型的响应进行仿真时，可使用可选的 x0 输入参量来指定非零初始状态值。假设有以下双状态 SISO 状态空间模型。

A = [-1.5 -3;
      3   -1];
B = [1.3; 0];
C = [1.15 2.3];
D = 0;
          
sys = ss(A,B,C,D);

假设您希望该系统在无输入的情况下从一组已知的初始状态开始演化 2 秒，然后施加单位阶跃变化。指定初始状态值向量 x0，并创建输入向量。

x0 = [-0.2 0.3];
t = 0:0.05:8;
u = zeros(length(t),1);
u(t>=2) = 1;
lsim(sys,u,t,x0)
grid on

该图的前半部分显示该系统从初始状态值 [-0.2 0.3] 开始的自由演化。在 t = 2 时，输入发生阶跃变化，该图显示该系统从该时间点的状态值开始对这个新信号的响应。

当您使用输出参量调用 lsim 时，它会以数组形式返回仿真响应数据。对于 SISO 系统，响应数据以列向量形式返回，该列向量与 t 长度相同。例如，提取 SISO 系统对方波的响应。使用 gensig 创建该方波。

sys = tf([2 5 1],[1 2 3]);
[u,t] = gensig("square",4,10,0.05);
[y,t] = lsim(sys,u,t);
size(y)

ans = 1×2

   201     1

向量 y 包含 t 中每个时间步时的仿真响应。（为方便起见，lsim 返回时间向量 t。）

对于 MIMO 系统，响应数据以 N×Ny×Nu 维度的数组形式返回，其中 Ny 和 Nu 是动态系统的输出数量和输入数量。例如，假设有以下状态空间模型，它表示具有两个输入和三个输出的三状态系统。

A = [-1.5  -0.2   1.0;
     -0.2  -1.7   0.6;
      1.0   0.6  -1.4];
  
B = [ 1.5  0.6;
     -1.8  1.0;
      0    0  ];

C = [ 0   -0.1 -0.2;
      0.7 -0.2 -0.3
     -0.65 0   -0.6];
  
D = [ 0.1  0;
      0.1  1.5
      0    0];

sys = ss(A,B,C,D);

提取三个输出通道对施加在两个输入端的该方波的响应。

uM = [u u];
[y,t] = lsim(sys,uM,t);
size(y)

ans = 1×2

   201     3

y(:,j) 是包含第 j 个输出对施加在两个输入端的该方波的响应的列向量。也就是说，y(i,:) 是一个包含三个值的向量，即第 i 个时间步的输出值。

由于 sys 是状态空间模型，因此您可以提取状态值随时间的演化，即状态值对输入信号的响应。

[y,t,x] = lsim(sys,uM,t);
size(x)

ans = 1×2

   201     3

x 的每一行包含在 t 中相应时间点的状态值 [x1,x2,x3]。换句话说，x(i,:) 是第 i 个时间步的状态向量。绘制状态值。

plot(t,x)

“绘制多个系统对同一输入的响应”示例说明了如何在一个轴上绘制多个单独系统的响应。当有多个动态系统排列在模型数组中时，lsim 会同时绘制所有这些系统的响应。

创建模型数组。对于此示例，请使用具有不同固有频率的二阶传递函数的一维数组。首先，为模型数组预分配内存。以下命令会创建一个行向量，其中包含 5 个零增益 SISO 传递函数。前两个维度表示模型输出和输入。其余维度是数组维度。（有关模型数组以及如何创建模型数组的详细信息，请参阅Model Arrays (Control System Toolbox)。）

sys = tf(zeros(1,1,1,5));

填充该数组。

w0 = 1.5:1:5.5;    % natural frequencies
zeta = 0.5;        % damping constant
for i = 1:length(w0)
   sys(:,:,1,i) = tf(w0(i)^2,[1 2*zeta*w0(i) w0(i)^2]);
end

绘制数组中所有模型对方波输入的响应。

[u,t] = gensig("square",5,15);
lsim(sys,u,t)

lsim 对数组中所有条目的响应使用相同的线型。区分各条目的一种方法是使用动态系统模型的 SamplingGrid 属性将数组中的每个条目与对应的 w0 值相关联。

sys.SamplingGrid = struct('frequency',w0);

现在，当您在 MATLAB® 图窗窗口中绘制响应时，您可以点击一条轨迹来查看它对应的频率值。

加载估计数据以估计模型。

load dcmotordata
z = iddata(y,u,0.1,'Name','DC-motor');

z 是一个 iddata 对象，它存储采样时间为 0.1 秒的单输入双输出估计数据。

使用估计数据 z 估计 4 阶状态空间模型。

[sys,x0] = n4sid(z,4);

sys 是估计得到的模型，x0 是估计得到的初始状态。

使用用于估计的相同输入数据和估计命令返回的初始状态，对 sys 的响应进行仿真。

[y,t,x] = lsim(sys,z.InputData,[],x0);

在本语法中，y 是系统响应，t 是用于仿真的时间向量，x 是状态轨迹。

比较两个输出的仿真响应 y 和测量响应 z.OutputData。

subplot(211), plot(t,z.OutputData(:,1),'k',t,y(:,1),'r')
legend('Measured','Simulated')
subplot(212), plot(t,z.OutputData(:,2),'k',t,y(:,2),'r')
legend('Measured','Simulated')

对于此示例，fcnMaglev.m 定义了一个磁悬浮系统的矩阵和偏移量。该磁悬浮系统通过线圈电流控制悬浮球的高度，线圈电流会对悬浮球产生磁力。此示例对模型进行开环仿真。

创建 LPV 模型。

lpvSys = lpvss('h',@fcnMaglev)

Continuous-time state-space LPV model with 1 outputs, 1 inputs, 2 states, and 1 parameters.
Model Properties

您可以使用圆点表示法设置额外的属性。

lpvSys.StateName = {'h','hdot'};
lpvSys.InputName = 'current';
lpvSys.InputName = 'height';

仿真此模型对任意正弦输入电流的响应。

h0 = 1;
[~,~,~,~,~,~,x0,u0,~] = fcnMaglev([],h0);
ic = findop(lpvSys,0,h0,x=x0,u=u0);
t = 0:1e-3:1;
u = u0*(1+0.1*sin(10*t));
y = lsim(lpvSys,u,t,ic,@(t,x,u) x(1));

绘制响应。

plot(t,y,t,u/u0)
legend('height','current')

当电流首次增大（h 减小）时，悬浮球会被磁铁吸引。随后电流的减小不足以使悬浮球回到原来的位置。

对于此模型，h = 0 是一个奇异点，即悬浮球碰到了磁铁。此时，LPV 模型不再有效。

采样时间的选择会极大地影响仿真结果。为了说明原因，假设有以下二阶模型。

w2 = 62.83^2;
sys = tf(w2,[1 2 w2]);

tau = 1;
Tf = 5;
Ts = 0.1;
[u,t] = gensig("square",tau,Tf,Ts);
lsim(sys,u,t)

figure
Ts2 = 0.01;
[u2,t2] = gensig("square",tau,Tf,Ts2);
lsim(sys,u2,t2)

此示例说明如何使用自定义仿真来获取状态快照数据并执行 POD 模型降阶。默认情况下，POD 算法提供三种类型的内置激励信号（啁啾信号、冲激信号和 PRBS 信号）来进行仿真。软件会对模型进行仿真并提取状态快照数据，然后近似计算可控性和可观测性 Gramian。或者，您可以提供通过 incrementalPOD 和 lsim 进行仿真生成的自定义 POD 数据。

生成具有 30 个状态、一个输入和一个输出的随机状态空间模型，并创建模型降阶任务。

rng(0)
sys = rss(30,1,1);
R = reducespec(sys,"pod");

对于此示例，创建一个叠加正弦信号作为输入信号来运行仿真。

t = linspace(0,100,10000);
u = 0.5*(sin(1.*t)+sin(3.*t)+sin(5.*t)+sin(8.*t)+sin(10.*t));

创建增量 POD 对象，用于存储可控性和可观测性 Gramian 的近似值。

rPOD = incrementalPOD;
oPOD = incrementalPOD;

使用自定义输入信号 u 对被控对象模型及其伴随模型进行仿真。

[~,~,~,~,rPOD] = lsim(sys,u,t,rPOD);
asys = adjoint(sys);
[~,~,~,~,oPOD] = lsim(asys,u,t,oPOD);

lsim 生成状态快照数据，并返回自定义 POD 数据作为输出。这些数据以与 R.Options 兼容的格式生成。

指定自定义数据并运行模型降阶算法。当您指定选项 CustomLr 和 CustomLo 时，软件会绕过内置仿真，直接使用这些数据。

R.Options.CustomLr = rPOD;
R.Options.CustomLo = oPOD;
R = process(R);

现在您便可以按照典型的工作流来选择阶数并获得降阶模型。

获得一个忽略总能量 0.01% 的降阶模型。

rsys = getrom(R,MaxLoss=1e-4);
order(rsys)

ans = 
9

bodeplot(sys,rsys,"r--")
legend("Original","Reduced")

ans = 
  Legend (Original, Reduced) with properties:

         String: {'Original'  'Reduced'}
       Location: 'northeast'
    Orientation: 'vertical'
       FontSize: 8.1000
       Position: [0.7707 0.8532 0.1875 0.0923]
          Units: 'normalized'

  Show all properties

降阶模型提供了全阶模型的良好逼近。

创建一个具有复系数的状态空间模型。

A = [-2-2i -2;1 0];
B = [2;0];
C = [0 0.5+2.5i];
D = 0;
sys = ss(A,B,C,D);

计算系统对周期为 4 秒的方波的响应。

[u,t] = gensig("square",4,12);
[y,tout] = lsim(sys,u,t);

得到的响应数据包含复数输出值。

y

y = 193×1 complex

   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
      ⋮