Matlab MPC toolbox tutorial 2

1.Introduction

MPC算法本身非常容易理解，关键在于如何将工程中的物理对象转换成合适的数学模型。

2.MPC Design

2.1时间参数

系统有控制周期$T_s$，prediction horizon $p$和control horizon $m$，这三个参数根据系统的开环属性进行设置，一般不会通过这三个参数调整系统特性。

• 控制周期
  控制频率应至少要大于控制对象带宽的4倍，控制周期越短越增加了控制器的运算量。理论上周期周期越短，对扰动的抑制作用越强，如果将控制周期缩短一倍，但是对相同的扰动的抑制，没有明显的增加，就不需要减短。
• prediction horizon
  如果期望的闭环响应时间是T，则T ≈ pTs
• control horizon
  一般control horizon 比较小1-3

2.2 constraints

添加hard constraints 可能会造成infeasible solution

• hard constraints
  MV 和MV Rate 的constraints 尽量不要全部设置成hard
  OV的constraints尽量全部设置成soft
• soft constraints ECR

0 — No violation allowed (hard constraint)
0.05 — Very small violation allowed (nearly hard)
0.2 — Small violation allowed (quite hard)
1 — average softness
5 — greater-than-average violation allowed (quite soft)
20 — large violation allowed (very soft)

• soft constraints weight
  比其他三项cost 之和要大1-2个数量级，如果violation 太大了，在这个基础上可以适当的扩大2-5倍。
• time variant constraints 和 time variant weight
  对于时变的限制和时变的权重，一般给prediction中的每个点分别分配限制和权重，但是这样调试起来及其麻烦。

2.3 调试权重

2.3.1调试工具：

• sensitivity:

[J,sens] = sensitivity(MPCobj,PerfFunc,PerfWeights,Tstop,r,v,simopt,utarget)

将主要项的cost记录下来，J是所有项的cost

• review
  将MPC控制做一个综述报告
  

2.3.2 weight tuning

$这个部分matlab还是有很多地方解释的非常不清楚：n是什么，nyc是什么$

2.4 plant with delay

• 设置模型
  如果模型满足下面的关系
  

% 用传递模型表示：只通过简单的输入输出的关系也可以作为系统模型
g11 = tf(12.8,[16.7 1],'IOdelay',1.0,'TimeUnit','minutes');
g12 = tf(-18.9,[21.0 1],'IOdelay',3.0,'TimeUnit','minutes');
g13 = tf(3.8,[14.9 1],'IOdelay',8.1,'TimeUnit','minutes');
g21 = tf(6.6,[10.9 1],'IOdelay',7.0,'TimeUnit','minutes');
g22 = tf(-19.4,[14.4 1],'IOdelay',3.0,'TimeUnit','minutes');
g23 = tf(4.9,[13.2 1],'IOdelay',3.4,'TimeUnit','minutes');
DC = [g11 g12 g13;
      g21 g22 g23];
DC.InputName = {'Reflux Rate','Steam Rate','Feed Rate'};
DC.OutputName = {'Distillate Purity','Bottoms Purity'};
DC = setmpcsignals(DC,'MD',3);
mpcDesigner(DC)

关于调试参数
$$P-M\gg t_{d,max}/\Delta t$$
$t_{d,max}$是最大延时时间；
$\Delta t$是控制周期

2.5 测试MPC控制器的鲁棒性

matlab的案例是比较，同一个控制器面对建模误差的鲁棒性。

• 模型误差

A = [-0.0285 -0.0014; -0.0371 -0.1476];
B = [-0.0850 0.0238; 0.0802 0.4462];
C = [0 1; 1 0];
D = zeros(2,2);
CSTR = ss(A,B,C,D);
CSTR.InputName = {'T_c','C_A_i'};
CSTR.OutputName = {'T','C_A'};
CSTR.StateName = {'C_A','T'};
CSTR = setmpcsignals(CSTR,'MV',1,'UD',2,'MO',1,'UO',2);

dB = [0 0;0.05 0];
perturbUp = CSTR;
perturbUp.B = perturbUp.B + dB;

perturbDown = CSTR;
perturbDown.B = perturbDown.B - dB;

• 观测step的反应

step(CSTR,perturbUp,perturbDown)
legend('CSTR','peturbUp','perturbDown')

曲线的差别：

• 在mpcDesigner中添加噪声，比较同一个控制器对不同模型的抗噪声的能力
• 

3.Applications

3.1 dc motor model

转换成物理模型：

matlab 上模型为：

模型状态变量是：

状态转移矩阵：

具体模型参数：输出参数$[\theta ,T]$，$\theta$是measured outputs, $T$是unmeasured outputs

最终的跟踪效果：

3.2 Paper Machine Process

化工造纸的过程：

状态变量x:


输入变量：$[G_p,G_w,N_p,N_w]$
控制变量$G_p,G_w$

Np(consistency of stock entering the feed tank)是measured disturbance,
Nw(the white water consistency)是unmeasured disturbance
输出变量: $[H_2,N_1,N_2]$
状态转移矩阵：

MPC模型：

跟踪的效果：

3.3倒立摆模型

倒立摆是典型的非线性模型：

• 线性模型的初始化：

opspec = operspec(mdlPlant);
opspec.States(1).Known = true;
opspec.States(1).x = 0;
opspec.States(3).Known = true;
opspec.States(3).x = 0;

• 将倒立摆模型进行线性化

io(1) = linio([mdlPlant '/dF'],1,'openinput');
io(2) = linio([mdlPlant '/F'],1,'openinput');
io(3) = linio([mdlPlant '/Pendulum and Cart System'],1,'openoutput');
io(4) = linio([mdlPlant '/Pendulum and Cart System'],3,'openoutput');
options = findopOptions('DisplayReport',false);
op = findop(mdlPlant,opspec,options);
plant = linearize(mdlPlant,op,io);
plant.InputName = {'dF';'F'};
plant.OutputName = {'x';'theta'};
pole(plant)
plant = setmpcsignals(plant,'ud',1,'mv',2);

mpc模型参数：