系列文章目录

---

前言

        本例演示如何通过命令行创建和测试模型预测控制器。

---

一、定义工厂模型

        本示例使用《使用 MPC Designer 设计控制器》中描述的工厂模型。创建工厂的状态空间模型，并设置一些可选的模型属性，如输入、状态和输出变量的名称和单位。

% continuous-time state-space matrices, with temperature as first output
A = [   -5  -0.3427;47.68    2.785];
B = [    0   10.3   0];
C = [0 1;1 0];
D = zeros(2,2);% create state space plant model
CSTR = ss(A,B,C,D);% set names
CSTR.InputName = {'T_c', 'C_A_f'};  % set names of input variables
CSTR.OutputName = {'T', 'C_A'};     % set names of output variables
CSTR.StateName = {'C_A', 'T'};      % set names of state variables% set units
CSTR.InputUnit = {'deg K', 'kmol/m^3'};     % set units of input variables
CSTR.OutputUnit = {'deg K', 'kmol/m^3'};    % set units of output variables
CSTR.StateUnit = {'kmol/m^3', 'deg K'};     % set units of state variables

        需要注意的是，该模型是由非线性模型围绕工作点线性化而来。因此，线性模型输入和输出信号的值代表与非线性模型中工作点值的偏差。更多信息，请参阅使用 MATLAB 代码进行线性化。

二、为不同的 MPC 类别分配输入和输出信号

        冷却剂温度是受控变量 (MV)，流入试剂浓度是未测量干扰输入 (UD)，反应器温度是测量输出 (MO)，试剂浓度是未测量输出 (UO)。

CSTR=setmpcsignals(CSTR,'MV',1,'UD',2,'MO',1,'UO',2);

三、显示工厂基本属性并绘制阶跃响应图

        使用阻尼来显示线性工厂模型极点的阻尼比、固有频率和时间常数。

damp(CSTR)

                                                                       Pole              Damping       Frequency      Time Constant  (rad/seconds)      (seconds)    -1.11e+00 + 1.09e+00i     7.13e-01       1.55e+00         9.03e-01    -1.11e+00 - 1.09e+00i     7.13e-01       1.55e+00         9.03e-01    

        绘制开环阶跃响应图。

step(CSTR)

        考虑到设备的额定稳定性，约 1 秒的时间常数建议采样时间不大于 0.5 秒。在采样时间为 0.5 秒的情况下，10 步的预测范围可以覆盖开环设备的整个稳定时间，因此可以使用这两个参数作为初始猜测。采样时间越短，意味着可用于控制计算的时间越短。预测时间越长（步数越多），意味着优化变量的数量越多，因此在可用时间步数内需要解决的计算问题也就越多。

四、创建控制器

        为提高示例的清晰度，请抑制来自 MPC 控制器的命令窗口消息。

old_status = mpcverbosity('off');

        创建一个模型预测控制器，控制间隔或采样时间为 0.5 秒，所有其他属性均为默认值，包括 10 步预测范围和 2 步控制范围。

mpcobj = mpc(CSTR,0.5) %#ok<*NOPTS>

MPC object (created on 12-Feb-2024 23:58:42):
---------------------------------------------
Sampling time:      0.5 (seconds)
Prediction Horizon: 10
Control Horizon:    2Plant Model:        --------------1  manipulated variable(s)   -->|  2 states  ||            |-->  1 measured output(s)0  measured disturbance(s)   -->|  2 inputs  ||            |-->  1 unmeasured output(s)1  unmeasured disturbance(s) -->|  2 outputs |--------------
Indices:(input vector)    Manipulated variables: [1 ]Unmeasured disturbances: [2 ](output vector)        Measured outputs: [1 ]Unmeasured outputs: [2 ]Disturbance and Noise Models:Output disturbance model: default (type "getoutdist(mpcobj)" for details)Input disturbance model: default (type "getindist(mpcobj)" for details)Measurement noise model: default (unity gain after scaling)Weights:ManipulatedVariables: 0ManipulatedVariablesRate: 0.1000OutputVariables: [1 0]ECR: 100000State Estimation:  Default Kalman Filter (type "getEstimator(mpcobj)" for details)UnconstrainedUse built-in "active-set" QP solver with MaxIterations of 120.

五、查看和修改控制器属性

        显示控制器属性及其当前值的列表。

get(mpcobj)

                          Ts: 0.5                 PredictionHorizon (P): 10                  ControlHorizon (C): 2                   Model: [1x1 struct]        ManipulatedVariables (MV): [1x1 struct]        OutputVariables (OV): [1x2 struct]        DisturbanceVariables (DV): [1x1 struct]        Weights (W): [1x1 struct]        Optimizer: [1x1 struct]        Notes: {}                  UserData: []                  History: 12-Feb-2024 23:58:42

        控制器显示的历史记录值会有所不同，因为这取决于控制器的创建时间。有关 MPC 控制器可编辑属性的说明，请在命令行输入 mpcprops。

        使用点符号修改这些属性。例如，将预测范围改为 15。

mpcobj.PredictionHorizon = 15;

        有些属性名称有别名。例如，可以使用别名 MV 代替 ManipulatedVariables。此外，许多控制器属性都是包含附加字段的结构。使用点符号可以查看和修改这些字段值。例如，由于默认情况下控制器对操纵变量和输出变量没有约束，因此可以使用点符号查看和修改这些约束。

        为控制器操作变量设置约束条件。

mpcobj.MV.Min = -10;    % K
mpcobj.MV.Max = 10;     % K
mpcobj.MV.RateMin = -1; % K/step
mpcobj.MV.RateMax = 1;  % K/step

        您可以缩写属性名称，但缩写必须明确无误。您还可以查看和修改控制器调整权重。例如，修改操作变量速率和输出变量的权重。

mpcobj.W.ManipulatedVariablesRate = 0.3;
mpcobj.W.OutputVariables = [1 0];

        您还可以在预测范围内定义随时间变化的约束条件和权重，预测范围在每个时间步长内 都会发生变化。例如，如果要在预测期结束时强制受控变量以更慢的速度变化，请输入

mpcobj.MV.RateMin = [-2; -1.5; -1; -1; -1; -0.5];mpcobj.MV.RateMax = [2; 1.5; 1; 1; 1; 0.5];

-0.5和 0.5 值用于第四步及以后。

        同样，您也可以为预测范围的每一步指定不同的输出变量权重。例如，输入

mpcobj.W.OutputVariables = [0.1 0; 0.2 0; 0.5 0; 1 0];

        您还可以修改控制器的干扰抑制特性。更多信息，请参阅 setEstimator、setindist 和 setoutdist。

六、审查控制器设计

        生成有关潜在运行时稳定性和性能问题的报告。

review(mpcobj)

        在本例中，审查命令发现设计存在两个潜在问题。第一个警告是由于 C_A 输出误差的权重为零。第二个警告是由于 MV 和 MVRate 都存在硬约束。

        您可以向下滚动查看每个测试结果的详细信息。

七、稳态闭环输出灵敏度增益

        计算闭环系统的闭环稳态增益矩阵。

SoDC = cloffset(mpcobj)

SoDC =-4.4409e-15

八、执行线性仿真

        使用 sim 函数对系统进行线性仿真。例如，仿真 mpcobj 在 26 个控制间隔内的闭环响应。从第二步开始，为反应器温度（第一输出）和试剂浓度（第二输出）分别指定 2 和 0 的设定点。请注意，由于第二输出的调整权重为 0，因此浓度的设定点将被忽略。

T = 26;
r = [0 0;2 0];
sim(mpcobj,T,r)

        您可以使用 mpcsimopt 修改仿真选项。例如，关闭操纵变量约束运行仿真。

mpcopts = mpcsimopt;
mpcopts.Constraints = 'off';
sim(mpcobj,T,r,mpcopts)

        现在，被操作变量的第一次移动超过了指定的 1 单位速率限制。

        您也可以在工厂/模型不匹配的情况下进行仿真。例如，定义一个工厂，其增益比控制器使用的模型中的增益大 50%，时间延迟为 0.1 秒。

mpcopts.Model = tf(1.5,1,'InputDelay',0.1)*CSTR;
sim(mpcobj,T,r,mpcopts)

        从闭环响应的振荡行为可以看出，工厂/模型不匹配会降低控制器的性能。退化可能很严重，必须逐个进行测试。

        其他仿真选项包括在被控变量或测量输出中添加指定的噪声序列、开环仿真和前瞻选项，以实现更好的设定点跟踪或测量干扰抑制。

九、存储和绘制仿真结果

        仿真系统，将仿真结果存储在 MATLAB® 工作区中。

[y,t,u] = sim(mpcobj,T,r);

        该语法抑制自动绘制，并以 y、t 和 u 变量形式返回仿真结果。您可以将结果用于其他目的，包括自定义绘图。例如，在同一图表中绘制操作变量和两个输出变量。

figuresubplot(2,1,1)
plot(t,u)
title('Inputs')
legend('T_c')subplot(2,1,2)
plot(t,y)
title('Outputs')
legend('T','C_A')
xlabel('Time')

        恢复 mpcverbosity 设置。

mpcverbosity(old_status);