预测控制matlab工具箱

1、matlab 预测控制工具箱一、设计示例注：示例均为matlab 自带1、伺服系统控制器的设计 图4-1 位置伺服系统原理图（1）伺服系统数学模型位置伺服系统由直流电机、变速箱、弹性轴、负载等组成（见图4-1）。可用如下微分方程来描述： 上式中的变量定义及取值见表4-1。欢迎流通！西安交通大学自动化系， 2将上述微分方程写成状态方程形式，有 其中，。（2）控制目标及约束控制目标：在电压V 的控制下，使过载角位置跟踪指定值。输出量仅有可观测。弹性轴承受的强度有一定的限制，因此对输出力矩T 的赋值作一定约束： 对输入电压的约束： 该系统有一个输入V ，两个输出：（可测量）和T （不可测量）。（3）

2、在matlab 中定义该系统的状态空间模型首先需要在matlab 中对系统的数学模型进行定义。可以直接在命令行输入mpcmotormodel （建议做法），也可以在命令行中输入下列命令：% DC-motor with elastic shaft%Parameters (MKS%-Lshaft=1.0; %Shaft lengthdshaft=0.02; %Shaft diametershaftrho=7850; %Shaft specific weight (Carbon steelG=81500*1e6; %Modulus of rigiditytauam=50*1e6; %Shear st

3、rengthMmotor=100; %Rotor massRmotor=.1; %Rotor radiusJmotor=.5*Mmotor*Rmotor2; %Rotor axial moment of inertia Bmotor=0.1; %Rotor viscous friction coefficient (A CASOR=20; %Resistance of armatureKt=10; %Motor constantgear=20; %Gear ratioJload=50*Jmotor; %Load inertiaBload=25; %Load viscous friction c

4、oefficientIp=pi/32*dshaft4; %Polar momentum of shaft(circular sectionKth=G*Ip/Lshaft; %Torsional rigidity(Torque/angleVshaft=pi*(dshaft2/4*Lshaft; %Shaft volumeMshaft=shaftrho*Vshaft; %Shaft massJshaft=Mshaft*.5*(dshaft2/4; %Shaft moment of inertiaJM=Jmotor;JL=Jload+Jshaft;Vmax=tauam*pi*dshaft3/16;

5、%Maximum admissible torqueVmin=-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*gear2-(Bmotor+Kt2/R/JM;3BB=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 modelsys=ss(AA,BB,Hyd;Hvd,Dy

6、d;Dvd;（4）利用MPCTOOL 界面设计控制器第一步：导入模型在命令行窗口中输入mpctool ，工具箱界面出来后点击Imort Plant，此时界面如4-2所示。 图4-2 Import Dialog Box with the Servomechanism Model Selected此时先前定义的伺服系统数学模型sys 默认被选中，Properties 列出了该模型的主要属性。点击Import ，导入数学模型，见4-3所示。第二步：设置属性系统的数学模型导入后，需要对输入变量和输出变量的属性做一些修改，使之符合实际情况。默认情况下，输入变量类型均是manipulated ，本例中不需

7、要修改。默认输出变量类型是measured ，本例的两个输出变量中，可测量，T 不可测量，需要相应修改。45 可以修改输入变量和输出变量的名称，增加描述性内容，设置单位和初始值。此处将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 interv

8、al(time units（采样周期）改为0.1，其它设置不变，见图4-4。 6图4-4 Controller Design View, Models and Horizons Pane接下来设置变量约束。点击Constraints 面板，输入变量V 、输出变量和T 的设置见图4-5。一般情况下需要对输入进行约束，而不对输出量作限制。角位置输出受到实际物理装置的限制，因此不作限制。对输出力矩T 设置了约束范围。 Max down rate 为非正值（或留空），它限制了每一个控制周期内控制量的减小幅度，同理，Max up rate为非负值（或留空），他对每一个控制周期内控制量的增幅做了限制。此处

9、均留空。 图4-5 Controller Design View, Constraints Pane 接下来点击Weight Tuning面板，此面板用于在设计控制器时对性能进行折中考虑。首先看Output weights，控制器的目的是为了尽可能缩小输出量与设定值之间的偏差。控制器计算预测时域内每个预测点处的实际值与设定值的误差，将之平方，然后相加，得到预测误差总和。控制器的一个目标就是使预测误差总和最小化，从而得到良好的跟踪效果。此处， 角位置输出应跟踪设定值，输出力矩T 则不需要有此限制。输入变量V 处于不断变化之中，Weight 应该设为0，为了使控制量变化平缓，此处设置Rate We

10、ight为0.1。 在进行权重设置时需要记住一点，权重间的相对大小的重要性要远远大于绝对大小。例如，变量A 变化0.1单位与变量B 变化100单位效果等同，那么在设置权重时，A 的权重将是B 权重的1000倍。 上方的滑动条用来整定控制器的总体性能。滑动条越靠近左方，对控制量的惩罚项越大，这会增加控制器的鲁棒性，但会降低跟踪精度。 当存在不可测量干扰时，Estimation 选项卡用来调节控制器的响应，此处没有用到。 7图4-6 Controller Design View, Weight Tuning Pane第三步：仿真场景的设置 见图4-7所示。Setpoints 用来设置输出参考信号。

11、设置ThetaL 的参考信号为阶跃信号，阶跃时间为1s ，幅度为1。设置T 的参考信号为常值信号，幅度为零。 Look ahead 选项假定参考信号的变化是已知的，这样控制器可以提前修正控制变量来提高跟踪精度，这在实际中不常见，此处不选中。 8图4-7 Simulation Settings View for Scenario1第四步：开始仿真 点击Simulate 按钮，或者工具栏中的绿色箭头，开始仿真。结果见4-8所示。灰色线为设定值，蓝色线为实际输出值。可见，输出信号在30s 时仍没有稳定下来，没有实现指令信号的准确跟踪。有一点值得称道，输出力矩和控制电压（此处没有给出）在指定范围之内。

12、0.511.52T h e t a L , R a d i a n s T , N . mPlant OutputsTime (sec图4-8 Response to Unit Step in the Angular Position Setpoint为了得到更快速、精确的响应，进入MPC1 Weight Tuning选项板，将滑动条拖至最右方，然后重新开始仿真，仿真结果见图4-9和图4-10。可见，此时的响应速度和跟踪精度都有了很大的改进。 00.811.21.4T h e t a L , R a d i a n s-80-60-40-20020T , N . mPlant OutputsT

13、ime (sec9图4-9 Faster Servo Response -5050100150200Plant Input: V, VTime (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 . mPlant OutputsTime (sec

14、图4-11 Servo Response for Step Increase of pi Radians10 -100-5050100150200250Plant Input: V, VTime (sec图4-12 Voltage Adjustments Corresponding to Figure 4-11第五步：保存工作 设计完毕后，将控制器导出至工作区，方便在Simulink block diagram进行分析和设计。 为了导出控制器，右击Controllers 节点，选择Export Controller，将会出现如图4-13所示的对话框，点击Export ，所设计的控制器将会以MP

15、C1为名导出到工作区。图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;OutputVa

16、riables(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 定义了权重值（

17、类似图4-6，此处定义的数值可以得到更好的性能）。Model 指定了待控制模型，Model.Plant.OutputGroup 指定第二个输出是不可测量的。输入下列语句，构建MPC 控制对象：ServoMPC = mpc(Model, Ts, PredictionHorizon, ControlHorizon;在构建MPC 控制对象后，可以对一些属性进行设置。例如可以利用下面语句之一对预测时域长度进行设置：ServoMPC.PredictionHorizon = 12;set(ServoMPC, PredictionHorizon, 12;为了得到对象的全部属性，输入： get(ServoMP

18、C为了得到对象的特定属性，例如控制时域长度，可以输入： M = get(ServoMPC, ControlHorizon;M = ServoMPC.ControlHorizon;可以同时设置多个属性。 本例中需要继续设置下列属性：set(ServoMPC, Weights, Weights, .ManipulatedVariables, ManipulatedVariables, .OutputVariables, OutputVariables;在命令行里直接输入属性名称可以显示属性值，也可以调用下列语句：display(ServoMPCsim命令用于线性系统的仿真：TimeSteps =

19、round(10/Ts;r = pi 0;y, t, u = sim(ServoMPC, TimeSteps, r;sim命令将输出结果和控制量分别用y 和u 来表示：subplot(311plot(t, y(:,1, 0 t(end, pi*1 1title(Angular Position (radians;subplot(312plot(t, y(:,2, 0 t(end, -78.5 -78.5title(Torque (nMsubplot(313stairs(t, utitle(Applied Voltage (voltsxlabel(Elapsed Time (seconds结果见

20、图4-14所示。Angular Position (radians Torque (nM Applied Voltage (voltsElapsed 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库，双击该模块，会出现如图

21、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 Scopes2、造纸机控制器的设计（1）造纸机非线性数学模型此例为matlab 预测控制工具箱自带，具体可参考Ying, Y., M.

22、 Rao, and Y. Sun Bilinear control strategy for paper making process, Chemical Engineering Communications (1992, Vol.111, pp. 13-28. 图4-18 Schematic of Paper Machine Headbox Elements设状态变量为，H 1 、H 2分别是Feed Tank和Headbox 液位高度，N 1 是feed tank浓度，N 2 是headbox 浓度。可测量输出为。 控制目标：使H 2和N 2跟踪设定值。 控制变量：。其中G p 是进入Fe

23、ed Tank液体流速，G w 是再回收水的流速。流入Feed Tank中液体浓度N p 为可测量干扰： 流入Feed Tank中再回收水的浓度N w 为不可测量干扰：变量均经过归一化处理。稳态时均为零，有界变量。时间单位为分。整个控制过程为开环稳定。mpcdemos 文件夹下包含了mpc_pmmodel.m文件，以S 函数形式定义了该控制过程的非线性数学模型。输入序列为G p , G w , N p , N w ，输出序列为H 2, N 1, N 2。（2）线性化原数学模型容易得到造纸机的head box线性化数学模型，得到的状态空间方程系数如下：A = -1.9300 0 0 00.394

24、0 -0.4260 0 00 0 -0.6300 00.8200 -0.7840 0.4130 -0.4260;B = 1.2740 1.2740 0 00 0 0 01.3400 -0.6500 0.2030 0.40600 0 0 0;C = 0 1.0000 0 00 0 1.0000 00 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

25、, N_2;（后两条指令用来在作图时标注图示，可选择）再利用下面语句获得该模型的阶跃响应，检查是否是开环稳定，结果见图4-19： step(PaperMach;-1-0.500.51From: Gp T o : H 2-1012 3T o : N 1 -10123T o : N2From: GwFrom: NpFrom: NwStep ResponseTime (secA m p l i t u d e图4-19 Linearized Paper Machine Models Step Responses从图4-19可以看出：（1）G p , G w 对三个输出均存在控制作用（存在耦合）。 （

26、2）G p , G w 对H 2的控制效果类似。 （3）N 2对G w 产生相反的响应过程。 这些特征使精确、独立地控制H 2 和N 2存在难度。（3）利用预测控制工具箱设计控制器在命令行中输入mpctool ，然后点击Import Plant，将先前创建的PaperMach 导入进来（具体过程可参考“伺服系统控制器的设计”第四小节）。 然后定义输入变量和输出变量的属性，最后得到的结果如图4-20所示。 图4-20 Signal Properties for the Paper Machine Application设置MPC 控制器参数（其它参数保留默认值）： Models and Hori

27、zons. 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. 对于 N 1, Weight = 0. (其它输出变量Weight = 1. 在Scenario1设置： Duration = 30 H 2 setpoint = 1 (constant 然后运行仿真，得到的结果见图4-21和

28、图4-22所示。00.511.5H 2, m-1.5-1-0.50N 1, % -0.10N 2, %Plant OutputsTime (sec图4-21 Servo Response for Unit Step in Headbox Level Setpoint -0.500.5G p , k g / h11.52G w , k g / h-101N p , %051015202530-101N w , %Plant InputsTime (sec图4-22 Manipulated Variable Moves Corresponding to Figure 4-21动态响应过程在8min

29、 左右，可以通过调整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峰值误差几乎降低了一个数量级，与此同

30、时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”选项框，然后运行，可以观察到在阶跃干扰未出现前，控制输入就存在

31、变化了，这是由于选中此栏后能够预知未来出现的干扰而提前作出变化。由于在实际中不可能提前得知未来出现的干扰，因此在下面的仿真中取消“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.6N 1,% -0.100.10.2N 2, %Plant OutputsTime (sec图4-23 Feedback Scenario: U