Matlab MPC toolbox tutorial 1

1.Introduction

MPC的基础原理并不复杂，在深入原理之前，借用Matlab MPC toolbox熟悉MPC控制架构。

1.1 MPC modeling

MPC关于模型的定义如下：

用simulink 表示上面的关系：

根据matlab对MPC处理的介绍：

The MPC controller performs all estimation and optimization calculations using a discrete-time, delay-free, state-space system with dimensionless input and output variables.

$MPC对输入输出的delay做了处理：$

Delay removal — If the discrete-time model includes any input, output, or internal delays, the absorbDelay command replaces them with the appropriate number of poles at z = 0, increasing the total number of discrete states. The InputDelay, OutputDelay, and InternalDelay properties of the resulting state-space model are all zero.

控制对象的模型（可以看成是开环模型）：

输入扰动模型（闭环模型的一部分）：

输出扰动模型（闭环模型的一部分）

测量噪声扰动模型（闭环模型的一部分）

1.2 MPC State Prediction

MPC是闭环控制模型，控制器状态，不仅只考虑开环模型的状态变量，还需要考虑各种扰动和测量误差，用于预测系统的下一个状态。

MPC控制器用于后续计算的状态方程稍不同于控制对象的状态方程：

• 控制器的状态变量由下面的状态变量组成：
  
• 控制器的输入变量：
  
• 更新后的状态矩阵如下：
  
• 加上状态观测器
  
  

1.3 quadratic program

MPC将线性模型转换成二次优化问题，基本的组成元素有：

cost的组成
在控制系统设计的时候，需要综合考虑各种指标，因而cost也是由多个项组合而成

• out reference tracking: 跟踪的精度
  
  
• Manipulated Variable Tracking
  输入权重，系统希望控制器输出的控制量越小越好（耗能更小…），如LQR控制器中的cost中有$uR{u}^T$项。下面的公式中，往往$u_{j,target}=0$。
  
  
• Manipulated Variable Move Suppression
  很多系统中倾向于控制变量的改变要比较smooth，尽量不出现较大的波动
  
  
• softening constraints
  $softening能解决特殊情况下需要突破constraints的情况$
  
  

完整的形式：

其中$V_{j,min}和V_{j,max}$是ECR value

QR的矩阵
假设状态转移矩阵如下：

转移矩阵如下：

在prediction horizon上y的预测数值：

简化这个公式：

系数矩阵为：

公式中${\Delta }_u$需要转换成控制变量$z_i$

举例：如果控制变量$z_i和{\Delta }_u$的关系如下图。

设置$\Delta u(0)=z0,\Delta u(2)=z1,\Delta u(5)=z2$

将状态方程进行简化：


化成二次优化的标准形式：

相关权重则表示为：

各种限制：

2. MPC Designer

2.1 流程

以matlab cstr(continuous stirred-tank reactor)[1]为例。反应外界环境的温度Tc是输入控制量，输入液体浓度$C_{Ai}$可能有一定的波动（并且没有办法去测量）；水箱中的温度T是可观测的输出量，$C_A$是无法观测的输出量；并且水箱温度T和$C_A$也是状态变量。


用一阶线性微分方程近似这个模型：

• step1:构建状态空间模型并设置模型具体结构
  
  设置具体结构：
  
• step2:定义每个状态的物理意义
  
• step3:设置限制和权重，hard constraints 和softening constraints 在二次优化方程中的形式如下。
  
  
• step4: 设置目标函数中的权重系数
  注意输出变量$C_A$因为观测变量，所以权重为0.
  
• step5: 设置输入信号，对控制器进行调试
  • 1.输入阶跃信号
    
  • 2.测试系统抗扰动的能力
    输入一个阶跃的输入扰动，系统很快适应下来。
    
• step6: 调试MPC控制器参数
  同样的根据[2],prediction horizon 定义process error的权重，control horizon 定义输入权重。
  本例中，control horizon 是 prediction horizon的$1/5$,在process weigh的设置的时候，给process error的权重0.3。
  

2.2 samples

2.2.1 单输入单输出样例

构建一个单输入单输出的二阶惯性模型：

在matlab里面，只要配置好相关的输入输出，模型和限制就可以。

plant = tf(1,[1 0 0]);
Ts = 0.1;   
p = 10;
m = 3;
mpcobj = mpc(plant, Ts, p, m);
mpcobj.MV = struct('Min',-1,'Max',5);   % 因为只有一个输入变量，只需要设置该变量的限制
mdl = 'mpc_doubleint';
open_system(mdl)
sim(mdl)

其他的内容：包括卡尔曼估计隐藏状态，以及matlab如何构建QP equation是看不到的。

对比PID控制器和MPC控制器的效果,单输入单输出模型，MPC的优势并不是太大。（下图中红色的是MPC，蓝色是PID）

比较MPC和PID控制器的输出指令，下图中（黄色的是PID 蓝色的是MPC）。

2.2.2 多输入单输出

假设多输入多输出的状态转移矩阵为：

% 设置MPC模型
sys = ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}));
Ts = 0.2;
model = c2d(sys,Ts);   
model = setmpcsignals(model,'MV',1,'MD',2,'UD',3);  % 3个输入变量，5个状态变量
% 设置MPC控制器
mpcobj = mpc(model,Ts,10,3);  % 控制器参数设置
mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);   % 设置输入限制
% 设置其他通道
mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);   % 噪声模型
% 设置sim
Tstop = 30;                                       % simulation time
Tf = round(Tstop/Ts);                        % number of simulation steps
r = ones(Tf,1);                                  % reference signal
v = [zeros(Tf/3,1);ones(2*Tf/3,1)];   % measured disturbance signal，如何作用在系统上
sim(mpcobj,Tf,r,v)
% 让仿真的效果更接近实际
SimOptions = mpcsimopt(mpcobj);   % 设置MPC的相关参数
d = [zeros(2*Tf/3,1);-0.5*ones(Tf/3,1)];        
SimOptions.Unmeas = d;                          % unmeasured input disturbance signal
SimOptions.OutputNoise=.001*(rand(Tf,1)-.5);    % output measurement noise
SimOptions.InputNoise=.05*(rand(Tf,1)-.5);      % noise on manipulated variables
% 手动修改kalman filter 用于估计无法直接观测的变量
[L,M,A1,Cm1] = getEstimator(mpcobj);
e = eig(A1-A1*M*Cm1)
poles = [.8 .75 .7 .85 .6 .81];
L = place(A1',Cm1',poles)';
M = A1\L;
setEstimator(mpcobj,L,M);

上面这个系统用simulink模拟如下图：

打印出模型的详细信息：

$输入扰动MD和UD是如何作用在这个系统上的$

2.2.3 非线性的多输入多输出问题

传递模型是一个非线性模型，采用EKF的思路，将当前时刻下的状态线性处理：

plant = linearize('mpc_nonlinmodel');
plant.InputName = {'Mass Flow';'Heat Flow';'Pressure'};  
plant.OutputName = {'Temperature';'Level'};
plant.InputUnit = {'kg/s' 'J/s' 'Pa'};
plant.OutputUnit = {'K' 'm'};
Ts = 0.2;                          
p = 5;
m = 2;
mpcobj = mpc(plant,Ts,p,m);
mpcobj.MV = struct('Min',{-3;-2;-2},'Max',{3;2;2},'RateMin',{-1000;-1000;-1000});
mpcobj.Weights = struct('MV',[0 0 0],'MVRate',[.1 .1 .1],'OV',[1 1]);
outdistmodel = tf({1 0;0 1},{[1 0 0 0],1;1,[1 0 0 0]});
setoutdist(mpcobj,'model',outdistmodel);
mdl2 = 'mpc_nonlinear_setoutdist';
open_system(mdl2)
sim(mdl2)

Reference

[1] https://www.youtube.com/watch?v=J6RqQ9qDDeU
[2] https://www.youtube.com/watch?v=gMOcBSmjdkQ