MATLAB控制系统建模（复习自用）

sys1=tf(num,den,Ts);

sys2=zpk(z,p,k,Ts);%零极点增益模型？

z-零点 p-极点 k-增益 Ts-采样时间

若无零点，则建立空矩阵占位 zpk([],[1,2],2,1)


[r,p,k]=residue(a,b);
[a,b]=residue(r,p,k);

状态空间ss

数学模型之间的转换

[A,B,C,D]=tf2ss(num,den) [z,p,k]=tf2zp(num,den)

[num,den]=ss2tf(A,B,C,D) [z,p,k]=ss2zp(A,B,C,D)

zp2tf zp2ss


连续系统到离散系统之间的转换

sysd=c2d(sys,Ts)%将连续系统转换为采样周期为Ts的离散系统

二、模型属性的定义

三、一个复杂系统的总的模型如何求取

% feedback parallel并联 series串联

% apeend connect

现在有如下图所示的模型，请用两种方法对其进行建模

clear,clc

num1=[1];den1=[10 1];
sys1=tf(num1,den1);
num2=[1];den2=[1 1];
sys2=tf(num2,den2);
num3=[1];den3=[1 2 1];
sys3=tf(num3,den3);
num4=[1];den4=[0.5 1];
sys4=tf(num4,den4);
num5=[3];den5=[4 1];
sys5=tf(num5,den5);
sys6=tf(1,1);%假想系统最后有一个增益为1 的环节

%% 方法一
s1=parallel(sys3,sys4);
s2=series(sys2,s1);
s3=feedback(s2,sys5);
s=series(sys1,s3)%TF模型

ss=ss(s)%转换为状态空间方程

%% 方法二
input=1;%选择系统1作为外部输入
output=6;%选择系统6作为外部输出
sys=append(sys1,sys2,sys3,sys4,sys5,sys6);
Q=[1 0 0;2 1 -5;3 2 0;4 2 0;5 3 4;6 3 4];%2 1 5表示第2个系统的输入为系统1和系统5
sysc=connect(sys,Q,input,output);%状态空间方程

sys=tf(sysc)%转换为TF模型

四、自己创建一个函数命令实现特定的功能

function的使用

例如：

1.创建一个函数，判断系统是否稳定

%% 该函数用于实现简单的判断系统是否稳定

function []=analysis_f(num0,den0)
% num0=[1 3];den0=[2 4 5 8 10];
% G=tf(num0,den0);
% Gc=feedback(G,1);
% [num1,den1]=tfdata(Gc)%返回值为元组cell
Gc=tf(num0,den0);
[num,den]=tfdata(Gc,'v');%返回值为向量，而非元组cell        !!!!!!!!!!!!!否则下列语句允许会出错！需要的是向量！
r=roots(den);
disp('系统的极点：');
disp(r);
% real(r);
a=find(real(r)>0);%寻找根在右半平面所在的位置%a是向量
% r(a) %寻找在右半平面的根
b=length(a);
if b>0
    disp('系统不稳定');
    fprintf('系统不稳定根的数目：%.d \n',b);
    disp('系统不稳定的根有:');% %.f',r(a));
    disp(r(a))
else
    disp('系统稳定')
end

2.创建一函数命令，分析系统的瞬时时域响应

%% 计算性能指标

function perfo(sys)
[y,t]=step(sys);
y_inf=y(end);%假设进入稳定状态时，系统最后一个值作为终值
ess=1-y_inf;

[y_max,max_idx]=max(y);
if y_max>y_inf*1.0001%有超调
    Mp=(y_max-y_inf)/y_inf*100;
    tp=t(max_idx);
    idx100=find(y>=y_inf,1);
    tr=t(idx100);
    tr_type='0% to 100%';
else%无超调
    Mp=0;
    tp=NaN;
    idx10=find(y>=0.1*y_inf,1);
    idx90=find(y>=0.9*y_inf,1);
    tr=t(idx90-idx10);
    tr_type='10% to 90%';
end

ts_idx=find(abs(y-y_inf)>0.02*y_inf,1,'last');
if ~isempty(ts_idx)
    ts=t(ts_idx+1);
else
    ts=t(1);
end

fprintf('--- 性能指标报告 ---\n');
fprintf('超调量 Mp: %.2f%%\n', Mp);
fprintf('上升时间 tr (%s): %.3f s\n', tr_type, tr);
fprintf('调节时间 ts (2%%): %.3f s\n', ts);
fprintf('稳态误差 ess: %.3f\n', ess);

end

五、如何使用S函数创建连续、离散或混合系统的仿真模型（S函数建模）

1.S函数的编写

①先说第一种：（我们老师教的）

function [sys,x0]=函数名(t,x,u,flag)
  if flag==1
      sys=连续系统一阶微分方程的值;
  elseif flag==2 %% 离散系统才有
      sys=离散系统下一刻的值;
  elseif flag==3
      sys=系统的输出值;
  elseif flag==0 %% 系统的初始设置参数
      sys=[连续状态的个数;离散系统状态的个数;系统输出个数;系统输入个数;...,
          系统的输出于输入有无直接关系;采样周期的个数];
      x0=状态初始值;%% e.g. x0=[0 0]
  else 
      sys[]=;
  end
end

②第二种是直接调用MATLAB内置的S函数模板

具体，调用如下：

在命令行窗口中输入

edit sfundsc2 %% 适用于2018版本及以上
%% 如若是2018版本以下的，听说是↓
edit sfuntempl %%末尾是l还是1我没试过，可以都试试（因版本原因，本人未操作过）

S函数模板：👇

function [sys,x0,str,ts,simStateCompliance] = sfundsc2(t,x,u,flag)
%SFUNDSC2 Example unit delay MATLAB File S-function
%   The MATLAB file S-function is an example of how to implement a unit
%   delay.
    
%   Copyright 1990-2009 The MathWorks, Inc.

switch flag,
 
  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,                                                
    [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes;    
    
  %%%%%%%%%%  
  % Update %
  %%%%%%%%%%
  case 2,                                               
    sys = mdlUpdate(t,x,u);
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                               
    sys = mdlOutputs(t,x,u);    

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case 9,                                               
    sys = [];

  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end

%end sfundsc2

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates  = 0;
sizes.NumDiscStates  = 1;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0  = 0;
str = [];
ts  = [.1 0]; % Sample period of 0.1 seconds (10Hz)

% speicfy that the simState for this s-function is same as the default
simStateCompliance = 'DefaultSimState';

% end mdlInitializeSizes

%
%=======================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=======================================================================
%
function sys = mdlUpdate(t,x,u)
sys = u;    

%end mdlUpdate

%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u)
sys = x;

%end mdlOutputs

将该.m文件的内容复制粘贴到你新创立的.m文件中，千万不要在原来的文件中更改！！！

该模板需要修改的地方有：

Ⅰ函数名

将sfunsc2改成你自己的函数名，注意要与你自己创立的.m文件名字要一样

Ⅱ修改初始化参数

在mdlInitializeSizes函数中，后面会举个例子

Ⅲ输出mdlOutputs

Ⅳ离散系统和连续系统有所不同

如若是连续系统，则没有case2，即没有更新参数这个情况，请删除case2这一分支

而若是离散系统，则需要case2

以下例子，是以连续系统进行说明的

Ⅴx1和x(1)是不同的

2.案例

3.考虑一个天线臂仰角控制系统，其数学模型为：

请利用DEE建模的方法、S函数的编写以及Simulink建模的方法给出当输入u为单位阶跃时x1和x2的变化波形。x1和x2的初始状态均为0（1，2请看成下标）。

这里只给出S函数的编写：

function [sys,x0,str,ts,simStateCompliance] = my_sfunction(t,x,u,flag)
%SFUNDSC2 Example unit delay MATLAB File S-function
%   The MATLAB file S-function is an example of how to implement a unit
%   delay.
    
%   Copyright 1990-2009 The MathWorks, Inc.

switch flag,
 
  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,                                                
    [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes;    
    
  %%%%%%%%%%  
  % Update %
  %%%%%%%%%%
  case 1,                                               
    sys = mdlUpdate(t,x,u);
    
  %%%%%%%%%%
  % Output %
  %%%%%%%%%%
  case 3,                                               
    sys = mdlOutputs(t,x,u);    

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case {2,9},                                               
    sys = [];

  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end
end
%end sfundsc2

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates  = 2;% 连续状态数：2个（x1=仰角，x2=角速度）
sizes.NumDiscStates  = 0;% 离散状态数：0（连续系统）
sizes.NumOutputs     = 2;% 输出数：2个（同时输出x1和x2）
sizes.NumInputs      = 1;% 输入数：1个（u=单位阶跃）
sizes.DirFeedthrough = 0;% 无直接馈通（输出不依赖输入u）
sizes.NumSampleTimes = 1;% 采样时间数量：1个

sys = simsizes(sizes);% 应用尺寸设置

x0  = [0,0];% 初始状态（实验要求：x1(0)=0，x2(0)=0）
str = [];% 预留参数
ts  = [0 0]; % Sample period of 0.1 seconds (10Hz)% 连续系统采样时间：[0 0]（无固定步长

% speicfy that the simState for this s-function is same as the default
simStateCompliance = 'DefaultSimState';

end
% end mdlInitializeSizes

%
%=======================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=======================================================================
%
function sys = mdlUpdate(t,x,u)
    x1 = x(1);  % 提取仰角
    x2 = x(2);  % 提取角速度

    % 代入非线性状态方程
    dx1dt = x2;                                 % dx1/dt = x2
    dx2dt = 9.81 * sin(x1) - 2 * x2 + u;        % dx2/dt = 9.81sin(x1) -2x2 + u

    sys = [dx1dt; dx2dt];  % 返回导数向量    
end
%end mdlUpdate

%
%=======================================================================
% mdlOutputs
% Return the output vector for the S-function
%=======================================================================
%
function sys = mdlOutputs(t,x,u)
sys = x;
end
%end mdlOutputs

S函数+simulink:

仿真结束后波形如下：

六、DEE(微分方程编辑器)工具箱对微分方程描述的系统进行建模

DEE工具箱当初我是费劲千辛万苦才在MATLAB社区上找到了它的相关包...

（等我弄明白怎么上传压缩包时再在此备份一份，工具箱导入的命令可以自行去问问AI）

然后我感觉我们老师好像也没怎么讲明白具体该怎么填写，也可能是我忘了，在此向老师致歉.......

还是上面五天线臂仰角控制系统的例子

如果你拥有了一个DEE工具箱，那么，请在你的命令行窗口输入：dee

（不区分大小写）并回车，你会得到一个类似于simulink的界面。聪明的你肯定已经学会了简单的模型建立↓（语句如有冒犯，真的很抱歉！！！）

双击DEE模块（即蓝色选中框），在页面中输入以下内容

仿真结果同使用S函数一样，就不在此浪费篇章了

七、利用SIMULINK建模

1.建模

欸，还是上面五天线臂仰角控制系统的例子

建立这种模型其实有个小技巧，就是你先确定变量有哪些，

像这种状态空间方程的式子，你需要先用积分模块integrator来建立导数与其自变量的关系

（二编：学了现代控制理论后，回过头来看，其实这就是状态空间表达式与系统结构模拟图的转换）

话不多说，先给出几个关键步骤：

双击simulink空白工作区，输入input+回车，output+回车，integrator+回车，建立上面的模型

由数学模型已知：

dx1=x2;

dx2=9.81*sin(x(1))-2*x(2)+u.

u这里题目说了是单位输入step

3.然后经过一些......就变成了下图所示的样子

2.To workspace

1.to workspace模块的用法

双击空白区域，输入 to workspace创立模块

①选择导出数据为时间序列

②数组

③结构体（无时间序列）

④结构体（带时间序列）

2.数据绘图

①时间序列

第一列是时间（out.x1.time），第二列是数据（out.x1.data）

% 1.时间序列
plot(out.x1.time,out.x1.data,out.x2.time,out.x2.data);
grid on

②数组

x1的数据👇

仿真时间👇

%% 2.数组
plot(out.tout,out.x1,out.tout,out.x2);
grid on;

③结构体（无时间序列，time为空矩阵）

%% 3.结构体（无时间序列）
plot(out.x1.signals.values,out.x2.signals.values);
grid on;
figure;
plot(out.tout,out.x1.signals.values,out.tout,out.x2.signals.values);
grid on;

八、ltiview 线性系统分析器

ltiview(sys);
ltiview(sys1,sys2);

方便修改图形

现在，要求使用4种不同方法给出系统的单位阶跃响应曲线

①lsim %%任意输入下的响应注意维数的匹配

②step %%直接绘图或用plot

③simulink

④ltiview

九、系统的优化设计

工具箱：sisotool (不做介绍)

在命令行窗口输入：sisotool

即可打开

系统的优化设计：寻找一个控制器参数使得目标函数的取值最优（最优指的是符合目标）

这里主要介绍单纯形寻优法

1.fminsearch %得出来的是变量，而非最小值

%使用无导数法计算无约束的多变量函数的最小值的非线性规划求解器。%%但容易只找到局部最优解

X=fminsearch(@sin,3)         %% 得出来的是变量，而非最小值 %3是初值，该函数依赖于初值的选取
X=fminsearch(@sin,1)

例题：

题外话：如果是ITAE则是这样↓（蓝色框框内的文字忘记修改了，其实这应该是ITAE准则）

IAE ISE ITAE ITSE ISTAE ISTSE准则 %% I积分，A绝对值，S平方，E误差, T时间

(从后往前看，可得到↓)

所以，IAE指的是对误差绝对值的积分；∫ || e || dt

ITSE指的是对误差的平方乘以时间后的积分 ∫ e^2*t dt

......

准则	全称	数学表达式	特点与适用场景
IAE	绝对误差积分	∫∥e(t)∥dt∫∥e(t)∥dt	最常用、最直观，对大小误差同等权重，易于实现。
ISE	平方误差积分	∫e(t)2dt∫e(t)2dt	更注重大误差（因为平方会放大），常用于强调需要快速抑制大偏差的场合（如航天、机器人），但可能导致响应振荡。
ITAE	时间乘以绝对误差积分	∫t∥e(t)∥dt∫t∥e(t)∥dt	惩罚持续的误差（时间t作为权重）。后期误差对指标影响更大，因此能有效减少调节时间，抑制超调后的长期波动，得到更平稳的响应。
ITSE	时间乘以平方误差积分	∫t⋅e(t)2dt∫t⋅e(t)2dt	ISE和ITAE的结合，同时强调大误差和持续误差，但计算和调参更复杂。

代码：

function ss = optm(x)
    global Kp Ki Kd;
    Kp = x(1); Ki = x(2); Kd = x(3);
    
    % 运行仿真
    yy = sim('PID_optimization.slx', 50);
    
    % 检查 yout 是否为空，防止报错
    if yy.yout.numElements > 0
        % 提取第一个输出端口（IAE 准则）的数据
        iae_values = yy.yout{1}.Values.Data;
        ss = iae_values(end); % 取最后一个积分值作为目标函数
    else
        error('Simulink 模型中未找到 Outport 模块，请检查 IAE 准则输出。');
    end
    
    % 容错处理：如果系统不稳定产生 NaN，给一个极大的惩罚值
    if isnan(ss) || isinf(ss)
        ss = 1e10; 
    end
end


%% 考试上机时使用下面代码↓（2018版本）

% [tt,xx,yy]=sim('PID_optimization.slx',50);%仿真时间50s

% ss=yy(end);

clear;
global Kp Ki Kd;
% 初始 PID 参数猜测 [Kp, Ki, Kd]
x0 = [1, 0.1, 0.01]; 

% 运行优化
% 'Display','iter' 可以让你在命令行实时看到误差 ss 正在下降
options = optimset('Display','iter','TolX',1e-3);
[best_x, min_iae] = fminsearch(@optm, x0, options);

fprintf('优化完成！\n最优参数：Kp=%.4f, Ki=%.4f, Kd=%.4f\n', best_x(1), best_x(2), best_x(3));
fprintf('迭代次数：%d \n',i);
fprintf('IAE准则下的最小值：%.4f \n',min_iae);

% 使用最优参数再跑一遍模型，观察对比效果
% 确保参数已更新
Kp = best_x(1); Ki = best_x(2); Kd = best_x(3);
out = sim('PID_optimization.slx', 50);

% 1. 绘制优化后的输出 (对应模型中的 out.y_optm 模块)
% 假设模块设置的是输出 timeseries 格式
plot(out.y_optm.Time, out.y_optm.Data, 'LineWidth', 1.5); 
hold on;

% 2. 绘制优化前的输出 (对应模型中的 out.y 模块)
plot(out.y.Time, out.y.Data, 'r--', 'LineWidth', 1); 
grid on;
legend('优化后 (y\_optm)', '优化前 (y)');
xlabel('时间 (s)');
ylabel('响应值');
title('PID 参数整定前后对比');

%% %% 实际考试在机房中用这个↓（2018版本）

%yy(end)

当前的 x 满足使用 1.000000e-03 的 OPTIONS.TolX 的终止条件，
F(X) 满足使用 1.000000e-04 的 OPTIONS.TolFun 的收敛条件

优化完成！
最优参数：Kp=12.9550, Ki=0.2036, Kd=15.9381
迭代次数：415
IAE准则下的最小值：5.0840

下面这个是之前简洁版本的（我们老师上课演示的代码），可以不用理会

结果：

1.初值为[100 0.1 1]：

Kp = 8.3013;

Ki = 0.1657;

Kd = -1.2509;

迭代次数（寻优次数）i=197;

根据IAE准则的最小值是6.6249