matlab练习程序（连续系统辨识与响应）

在不知道系统模型的情况下，如果想要使用模型的方法控制系统，可以先用系统辨识求得系统的状态空间参数，然后设计控制器。

下面根据已知状态参数设计了一个系统，再利用控制量得到系统响应。

然后结合控制量和系统响应，利用matlab自带的子空间迭代方法辨识得到系统的状态方程。

再比较一下原始响应和辨识结果响应。

又比较了不同方法下的系统的冲击与阶跃响应。

matlab代码如下：

clear all;close all;clc;

load dryer2
warning off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%连续系统的系统辨识%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = [-1.5,-2;1,0];
B = [0.5;0];
C = [0,1];
D = 0;
sys = ss(A,B,C,D);              %状态方程
sys_tf = tf(sys);               %传递函数

u = u2(1:100);                  %控制量
y = lsim(sys,u,(1:100)');       %响应

idata = iddata(y,u);            %构造辨识结构体
idsys = n4sid(idata,2,'Ts',0,'Form','canonical'); %系统辨识得到状态方程(连续)
idsys_tf = idtf(idsys);         %构造新的传递函数

y2 = lsim(idsys,u,(1:100)');    %用原始控制量控制已辨识的系统生成响应

plot(y,'r-*');
hold on;
plot(y2,'go');
legend('原始响应','利用辨识结果得到的响应');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%状态空间转传递函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms s;
Gc = sys.C*inv(s*eye(2) - sys.A)*sys.B+sys.D
sys_tf

%%%%%%%%%%%%对传递函数计算反拉普拉斯变换并和原系统比较冲击和阶跃响应%%%%%%%%%%
y = ilaplace(Gc);
evaly = @(t) vpa(eval(y),10);
figure;
u=0:0.1:10;
imp = double(evaly(u));
subplot(4,2,1);
plot(u,imp,'r');
subplot(4,2,2);
ste = cumsum(imp)/10.0;
plot(u,ste,'b');
subplot(4,2,3);
[r,t] = impulse(sys);
plot(t,r,'r')
subplot(4,2,4);
[r,t] = step(sys);
plot(t,r,'b');

%%%%%%%%%%%%%%利用符号方法解上述系统的微分方程并求阶跃和冲击响应%%%%%%%%%%%%%%
y = dsolve('D2y+1.5*Dy+2*y=0.5','Dy(0)=0','y(0)=0','x');
evaly = @(x) vpa(eval(y),10);
x=0:0.1:10;
imp = double(evaly(x));
subplot(4,2,5);
plot(x(1:end-1),diff(imp)./diff(x),'r');
subplot(4,2,6);
plot(x,imp,'b');

%%%%%%%%%%%%%%%%%%利用数值方法解上述微分方程并求阶跃和冲击响应%%%%%%%%%%%%%%%%
[x,y] = ode45(@func,[0 10],[0; 0]);
subplot(4,2,7);
plot(x(1:end-1),diff(y(:,2))./diff(x),'r')
subplot(4,2,8);
plot(x,y(:,2),'b')

function dy = func(x,y)
dy = [y(2);
    0.5*x-1.5*y(2)-2*y(1)];
end

结果如下：

用原控制量与辨识结果生成响应：

不同方法下得到的冲击与阶跃响应：