预测控制matlab工具箱

matlab预测控制工具箱

一、设计示例

注：示例均为matlab自带

1、伺服系统控制器的设计

图4-1 位置伺服系统原理图

（1）伺服系统数学模型

位置伺服系统由直流电机、变速箱、弹性轴、负载等组成（见图4-1）。可用如下微分方程来描述：

上式中的变量定义及取值见表4-1。

将上述微分方程写成状态方程形式，有

其中，。

（2）控制目标及约束

控制目标：在电压V的控制下，使过载角位置跟踪指定值。输出量仅有可观测。

弹性轴承受的强度有一定的限制，因此对输出力矩T的赋值作一定约束：

对输入电压的约束：

该系统有一个输入V，两个输出：（可测量）和T（不可测量）。

（3）在matlab中定义该系统的状态空间模型

首先需要在matlab中对系统的数学模型进行定义。可以直接在命令行输入mpcmotormodel（建议做法），也可以在命令行中输入下列命令：

% DC-motor with elastic shaft

%

%Parameters (MKS)

%-----------------------------------------------------------

Lshaft=1.0; %Shaft length

dshaft=0.02; %Shaft diameter

shaftrho=7850; %Shaft specific weight (Carbon steel)

G=81500*1e6; %Modulus of rigidity

tauam=50*1e6; %Shear strength

Mmotor=100; %Rotor mass

Rmotor=.1; %Rotor radius

Jmotor=.5*Mmotor*Rmotor^2; %Rotor axial moment of inertia Bmotor=0.1; %Rotor viscous friction coefficient (A CASO)

R=20; %Resistance of armature

Kt=10; %Motor constant

gear=20; %Gear ratio

Jload=50*Jmotor; %Load inertia

Bload=25; %Load viscous friction coefficient

Ip=pi/32*dshaft^4; %Polar momentum of shaft

(circular) section

Kth=G*Ip/Lshaft; %Torsional rigidity

(Torque/angle)

Vshaft=pi*(dshaft^2)/4*Lshaft; %Shaft volume

Mshaft=shaftrho*Vshaft; %Shaft mass

Jshaft=Mshaft*.5*(dshaft^2/4); %Shaft moment of inertia

JM=Jmotor;

JL=Jload+Jshaft;

Vmax=tauam*pi*dshaft^3/16; %Maximum admissible torque

Vmin=-Vmax;

%Input/State/Output continuous time form

%----------------------------------------------------------

AA=[0 1 0 0;

-Kth/JL -Bload/JL Kth/(gear*JL) 0;

0 0 0 1;

Kth/(JM*gear) 0 -Kth/(JM*gear^2)

-(Bmotor+Kt^2/R)/JM];

BB=[0;0;0;Kt/(R*JM)];

Hyd=[1 0 0 0];

Hvd=[Kth 0 -Kth/gear 0];

Dyd=0;

Dvd=0;

% Define the LTI state-space model

sys=ss(AA,BB,[Hyd;Hvd],[Dyd;Dvd]);

（4）利用MPCTOOL界面设计控制器

第一步：导入模型

在命令行窗口中输入mpctool，工具箱界面出来后点击Imort Plant…，此时界面如4-2所示。

图4-2 Import Dialog Box with the Servomechanism Model Selected 此时先前定义的伺服系统数学模型sys默认被选中，Properties列出了该模型的主要属性。

点击Import，导入数学模型，见4-3所示。

第二步：设置属性

系统的数学模型导入后，需要对输入变量和输出变量的属性做一些修改，使之符合实际情况。默认情况下，输入变量类型均是manipulated，本例中不需要修改。默认输出变量类

型是measured，本例的两个输出变量中，可测量，T不可测量，需要相应修改。

可以修改输入变量和输出变量的名称，增加描述性内容，设置单位和初始值。此处将In1,Out1,Out2修改为V,ThetaL,T，其它修改见图4-3。

图4-3 Design Tool After Importing the Plant Model and Specifying Signal Properties 本例中的模型为连续系统，需要将之变为离散系统。选中MPCdesign---Controllers---MPC1---Model and Horizons，将Control interval(time units)（采样周期）改为0.1，其它设置不变，见图4-4。

图4-4 Controller Design View, Models and Horizons Pane 接下来设置变量约束。点击Constraints面板，输入变量V、输出变量和T的设置见图4-5。一般情况下需要对输入进行约束，而不对输出量作限制。角位置输出受到实际物

理装置的限制，因此不作限制。对输出力矩T设置了约束范围。

rate为非正值（或留空），它限制了每一个控制周期内控制量的减小幅度，同down

Max

理，Max up rate为非负值（或留空），他对每一个控制周期内控制量的增幅做了限制。此处均留空。

图4-5 Controller Design View, Constraints Pane

接下来点击Weight Tuning面板，此面板用于在设计控制器时对性能进行折中考虑。首先看Output weights，控制器的目的是为了尽可能缩小输出量与设定值之间的偏差。控制器计算预测时域内每个预测点处的实际值与设定值的误差，将之平方，然后相加，得到预测误差总和。控制器的一个目标就是使预测误差总和最小化，从而得到良好的跟踪效果。此处，

角位置输出应跟踪设定值，输出力矩T则不需要有此限制。

输入变量V处于不断变化之中，Weight应该设为0，为了使控制量变化平缓，此处设置Rate Weight为0.1。

在进行权重设置时需要记住一点，权重间的相对大小的重要性要远远大于绝对大小。例如，变量A变化0.1单位与变量B变化100单位效果等同，那么在设置权重时，A的权重将是B权重的1000倍。

上方的滑动条用来整定控制器的总体性能。滑动条越靠近左方，对控制量的惩罚项越大，这会增加控制器的鲁棒性，但会降低跟踪精度。

当存在不可测量干扰时，Estimation选项卡用来调节控制器的响应，此处没有用到。

图4-6 Controller Design View, Weight Tuning Pane

第三步：仿真场景的设置

见图4-7所示。Setpoints用来设置输出参考信号。设置ThetaL的参考信号为阶跃信号，阶跃时间为1s，幅度为1。设置T的参考信号为常值信号，幅度为零。

Look

ahead选项假定参考信号的变化是已知的，这样控制器可以提前修正控制变量来提

高跟踪精度，这在实际中不常见，此处不选中。

图4-7 Simulation Settings View for "Scenario1"

第四步：开始仿真 点击Simulate 按钮，或者工具栏中的绿色箭头，开始仿真。结果见4-8所示。灰色线为设定值，蓝色线为实际输出值。可见，输出信号在30s 时仍没有稳定下来，没有实现指令信号的准确跟踪。有一点值得称道，输出力矩和控制电压（此处没有给出）在指定范围之内。

0.511.52T h e t a L , R a d i a n

s

T , N .m

Plant Outputs

Time (sec)

图4-8 Response to Unit Step in the Angular Position Setpoint

为了得到更快速、精确的响应，进入MPC1 Weight Tuning 选项板，将滑动条拖至最右方，然后重新开始仿真，仿真结果见图4-9和图4-10。可见，此时的响应速度和跟踪精度都有了很大的改进。

00.20.40.60.811.21.4T h e t a L , R a d i a n

s

-80

-60-40-20020T , N .m

Plant Outputs

Time (sec)

图4-9 Faster Servo Response

-50

50

100

150

200

Plant Input: V, V

Time (sec)

图4-10 Manipulated Variable Adjustments Corresponding to Figure 4-9

将ThetaL 的阶跃响应幅度重新设置为pi ，仿真结果见图4-11和图4-12所示。可见仍能实现很好的跟踪效果，注意此时对输出力矩和控制电压做了限幅。

01234T h e t a L , R a d i a n s

-80

-60-40-2002040T , N .m

Plant Outputs

Time (sec)

图4-11 Servo Response for Step Increase of pi Radians

-100

-50

50

100

150

200

250

Plant Input: V, V

Time (sec)

图4-12 Voltage Adjustments Corresponding to Figure 4-11

第五步：保存工作 设计完毕后，将控制器导出至工作区，方便在Simulink block diagram 进行分析和设计。 为了导出控制器，右击Controllers 节点，选择Export Controller ，将会出现如图4-13所示的对话框，点击Export ，所设计的控制器将会以MPC1为名导出到工作区。

图4-13 Exporting a Controller to the Workspace

在退出MPCTOOL 时，会提示将整个设计保存到一个MAT 文件中。

（5）利用MPC Toolbox 命令函数设计控制器

熟悉MPC 工具箱后，也许你会觉得利用命令函数设计控制器更为方便。

在步骤（3）定义的数学模型的基础上，输入下列指令：

ManipulatedVariables = struct('Min', -220, 'Max', 220, 'Units', 'V');

OutputVariables(1) = struct('Min', -Inf, 'Max', Inf, 'Units',

'rad');

OutputVariables(2) = struct('Min', -78.5, 'Max', 78.5, 'Units',

'Nm');

Weights = struct('Input', 0, 'InputRate', 0.05, 'Output', [10 0]); Model.Plant = sys;

Model.Plant.OutputGroup = {[1], 'Measured' ; [2], 'Unmeasured'};

Ts = 0.1;

PredictionHorizon = 10;

ControlHorizon = 2;

这样就生成了一些结构变量。ManipulatedVariables定义了控制电压的变化范围、单位等，Weights定义了权重值（类似图4-6，此处定义的数值可以得到更好的性能）。Model指定了待控制模型，Model.Plant.OutputGroup指定第二个输出是不可测量的。

输入下列语句，构建MPC控制对象：

ServoMPC = mpc(Model, Ts, PredictionHorizon, ControlHorizon);

在构建MPC控制对象后，可以对一些属性进行设置。例如可以利用下面语句之一对预测时域长度进行设置：

ServoMPC.PredictionHorizon = 12;

set(ServoMPC, 'PredictionHorizon', 12);

为了得到对象的全部属性，输入：

get(ServoMPC)

为了得到对象的特定属性，例如控制时域长度，可以输入：

M = get(ServoMPC, 'ControlHorizon');

M = ServoMPC.ControlHorizon;

可以同时设置多个属性。

本例中需要继续设置下列属性：

set(ServoMPC, 'Weights', Weights, ...

'ManipulatedVariables', ManipulatedVariables, ...

'OutputVariables', OutputVariables);

在命令行里直接输入属性名称可以显示属性值，也可以调用下列语句：

display(ServoMPC)

sim命令用于线性系统的仿真：

TimeSteps = round(10/Ts);

r = [pi 0];

[y, t, u] = sim(ServoMPC, TimeSteps, r);

sim命令将输出结果和控制量分别用y和u来表示：

subplot(311)

plot(t, y(:,1), [0 t(end)], pi*[1 1])

title('Angular Position (radians)');

subplot(312)

plot(t, y(:,2), [0 t(end)], [-78.5 -78.5])

title('Torque (nM)')

subplot(313)

stairs(t, u)

title('Applied Voltage (volts)')

xlabel('Elapsed Time (seconds)')

结果见图4-14所示。

Angular Position (radians)

Torque (nM)

Applied Voltage (volts)

Elapsed Time (seconds)

图4-14 Plotting the Output of the sim Command

（6）在Simulink中使用MPC Tools

图4-15为在Simulink中搭建的伺服控制系统的预测控制器。大部分模块来自Simulink 常用库。伺服系统的数学模型来自Control System Toolbox库，sys必须在工作区间存在定义。MPC Controller来自MPC Blocks库，双击该模块，会出现如图4-16所示的对话框，在MPC Controller中输入ServoMPC（ServoMPC必须存在于工作区间，其创建过程可参考上一节）。

仿真结果见图4-17。

图4-15 Block Diagram for the Servomechanism Example

图4-16 Model Predictive Control Toolbox? Simulink? Block Dialog Box

图4-17 Servomechanism Simulation Scopes

2、造纸机控制器的设计

（1）造纸机非线性数学模型

此例为matlab预测控制工具箱自带，具体可参考Ying, Y., M. Rao, and Y. Sun "Bilinear control strategy for paper making process," Chemical Engineering Communications (1992), V ol. 111, pp. 13-28.

图4-18 Schematic of Paper Machine Headbox Elements

设状态变量为，H1、H2分别是Feed Tank和Headbox液位高度，N1是feed tank浓度，N2是headbox浓度。

可测量输出为。

控制目标：使H2和N2跟踪设定值。

控制变量：。其中G p是进入Feed Tank液体流速，G w是再回收水的流速。

流入Feed Tank中液体浓度N p为可测量干扰：

流入Feed Tank中再回收水的浓度N w为不可测量干扰：

变量均经过归一化处理。稳态时均为零，有界变量。时间单位为分。整个控制过程为开环稳定。

mpcdemos文件夹下包含了mpc_pmmodel.m文件，以S函数形式定义了该控制过程的非线性数学模型。输入序列为G p, G w, N p, N w，输出序列为H2, N1, N2。

（2）线性化原数学模型

容易得到造纸机的head box线性化数学模型，得到的状态空间方程系数如下：

A = [-1.9300 0 0 0

0.3940 -0.4260 0 0

0 0 -0.6300 0

0.8200 -0.7840 0.4130 -0.4260];

B = [1.2740 1.2740 0 0

0 0 0 0

1.3400 -0.6500 0.2030 0.4060

0 0 0 0];

C = [0 1.0000 0 0

0 0 1.0000 0

0 0 0 1.0000];

D = zeros(3,4);

然后利用下述命令创建一个连续时间LTI状态空间模型：

PaperMach = ss(A, B, C, D);

PaperMach.InputName = {'G_p', 'G_w', 'N_p', 'N_w'};

PaperMach.OutputName = {'H_2', 'N_1', 'N_2'};

（后两条指令用来在作图时标注图示，可选择）

再利用下面语句获得该模型的阶跃响应，检查是否是开环稳定，结果见图4-19：step(PaperMach);

-1-0.500.51From: G

p

T o : H 2

-1012

3T o : N 1

-10123T o : N

2

From: G

w

From: N

p

From: N

w

Step Response

Time (sec)

A m p l i t u d e

图4-19 Linearized Paper Machine Model's Step Responses

从图4-19可以看出：

（1）G p , G w 对三个输出均存在控制作用（存在耦合）。 （2）G p , G w 对H 2的控制效果类似。 （3）N 2对G w 产生相反的响应过程。 这些特征使精确、独立地控制H 2 和N 2存在难度。

（3）利用预测控制工具箱设计控制器

在命令行中输入mpctool ，然后点击Import Plant…，将先前创建的PaperMach 导入进来（具体过程可参考“伺服系统控制器的设计”第四小节）。 然后定义输入变量和输出变量的属性，最后得到的结果如图4-20所示。

图4-20Signal Properties for the Paper Machine Application 设置MPC控制器参数（其它参数保留默认值）：

?Models and Horizons. Control interval = 2 minutes

?Constraints. 对于G p和G w, Minimum = -10, Maximum = 10, Max down rate = -2, Max up rate = 2.

?Weight Tuning. 对于G p和G w, Weight = 0, Rate weight = 0.4.

对于N1, Weight = 0. (其它输出变量Weight = 1.)

在Scenario1设置：

? Duration = 30

?H2 setpoint = 1 (constant)

然后运行仿真，得到的结果见图4-21和图4-22所示。

00.511.5H 2,

m

-1.5-1-0.50N 1,

%

-0.1

00.10.20.3N 2, %

Plant Outputs

Time (sec)

图4-21 Servo Response for Unit Step in Headbox Level Setpoint

-0.500.5G p , k g /

h

11.52G w , k g /

h

-101N p ,

%

0510********

-1

01N w , %

Plant Inputs

Time (sec)

图4-22 Manipulated Variable Moves Corresponding to Figure 4-21

动态响应过程在8min 左右，可以通过调整control interval 、输入变量的rate weights 、up/down rate constraints 来减小动态响应时间。为了获得鲁棒性，前面控制器参数设置得较

保守。在现在控制器的基础上，继续改进。 由图4-21可见N 1存在稳态误差（-0.25单位）。由于只存在两个输入变量，不可能使三个输出变量同时跟踪设定值。此处不需要跟踪N 1，所以设置它的weight 为零，否则会使三个输出均存在稳态误差（可以将weight 改为1来验证）。 浓度控制的重要性要高于液面的控制。将H 2 weight 由1改为0.2，仿真后会发现N 2峰值误差几乎降低了一个数量级，与此同时H 2的响应时间由8min 延长到18min 。（在后面的仿真中继续使用此组参数） 点击Scenarios ，新建一个名为Feedforward 的场景。在新建场景中，令可测量干扰N p 作阶跃变化：Initial value = 0, Size = 1, Time = 10。令所有输出变量的参考值为零，仿真时间为30min 。运行后，会发现H 2 和N 2离设定值有轻微偏离。下面使用“look ahead ”功能。注意在先前的仿真中，阶跃干扰出现以前（10min ），控制输入变量一直为零。选中N p 栏中的“look ahead ”选项框，然后运行，可以观察到在阶跃干扰未出现前，控制输入就存在变化了，这是由于选中此栏后能够预知未来出现的干扰而提前作出变化。由于在实际中不可能提前得知未来出现的干扰，因此在下面的仿真中取消“look ahead ”选项。 新建一个名为“Feedback ”的仿真场景，观测对不可测量干扰N w 的响应性能。将N p 设为0（constant ），将N w 设为1（constant ），表示在时间0时突然出现一个幅值为1的未知干扰。运行后得到如图4-23的仿真结果。H 2 和N 2离设定值0有轻微偏差。与图4-21相比，响应时间较长。Model and Horizons 选项卡中的control interval 会影响控制性能。估计器的设置也会影响对未知干扰的控制性能。图4-23为Estimation 默认设置时的结果。

-0.2-0.15-0.1-0.050H 2,

m

-0.200.20.40.6N 1,

%

-0.1

00.10.2N 2, %

Plant Outputs

Time (sec)

图4-23 Feedback Scenario: Unmeasured Disturbance Rejection

选中MPC1---Estimation ，出现如图4-24所示的对话框。最底下的状态栏提示信息表示使用预测控制工具箱默认参数。默认进入Output Disturbances 选项卡，Signal-by-signal 被默认选中。前两个输出量设为独立阶跃干扰。然后点击Input Disturbances ，不可测输入干扰也为独立阶跃类型。这样，一共有三个独立阶跃类型的干扰，幅值均为1。选中

Measurement Noise选项卡，看到每个输出量均被加上白噪声。如果某个输出干扰很大，那么控制器就会减小该输出的权重，依靠该输出的预测值。在此处造纸机控制器的设计中，默认的噪声设置是合理的。

图4-24 Default Estimator Assumptions: Output Disturbances

（4）利用Simulink对原非线性模型进行仿真

上一节中利用线性化模型设计了预测控制器，得到了很好的控制效果。这样子得到的控制器没有预测误差。由于实际模型是非线性的，需要将得到的控制器控制非线性模型，这种情况下会引入预测误差。

图4-25为造纸机非线性数学模型的Simulink控制框图，图中的S函数为造纸机的非线性数学模型，其定义见mpc_pmmodel.m。