第7章模型预测控制系统的计算机辅助设计

1、第7章 模型预测控制系统的 计算机辅助设计 7.1 模型预测控制工具箱概述7.2 基于阶跃响应的模型预测控制7.3 基于状态空间模型的模型预测控制7.1 模型预测控制工具箱概述 模型预测控制（MPC）工具箱是分析和设计模型预测控制系统的函数集合。 模型预测控制首先于20世纪70年代运用到工业控制中, 到80年代后期得到广泛应用。 现在, 该方法已经广泛应用于化工工业和其它工业领域的多变量控制系统。 MPC方法几乎可以适用于任何控制问题, 它主要用来解决以下类型的问题： (1) 系统具有大量需要控制的变量。 (2) 被控变量具有约束。 (3) 控制目标的改变以及控制/传感设备的失效。 (4) 具

2、有时延的系统。 在实际工程中, MPC方法又经常被分为动态矩阵控制（DMC）、 IDCOM和模型算法控制等等。 7.2 基于阶跃响应的模型预测控制 7.2.1 阶跃响应模型 对于一个线性时不变（LTI）的单输入单输出（SISO）系统, 设系统在单位输入变化v下的系统输出变化为 0, , s1, s2, , sn, 这里我们假设系统在n步以后进入稳定状态, 这样, 系统开始n步s1, s2, , sn的阶跃响应可以完全刻画系统的模型, 我们可以利用这些阶跃响应计算任何输入序列下的系统输出响应： (7.1) 上述阶跃响应模型也可以用于稳定且具有积分性质的过程。 对于一个积分过程, 我们可以假设系统

3、在n步以后的响应信号的斜率保持常数, 也就是说： (7.2) 对于具有nv个输入和ny个输出的多输入多输出（MIMO）系统, 我们可以得到一系列阶跃响应的系数矩阵(7.3) 其中, sl,m,i是第m个输入到第l个输出的第i个阶跃响应。 MPC工具箱将按照下面的格式存储阶跃响应模型： 其中的delt2是系统采样时间; 向量nout表示相应的输出是否正在进行积分： nout(i) = 0表示第i个输出正在积分， nout(i) = 1表示第i个输出是稳定的。 系统的阶跃响应可以直接从系统辨识数据中获得, 也可以由连续或离散传递函数和状态空间模型产生。 例如, 如果某个离散系统的描述（采样时间假定

4、为T = 0.1 s）为 下面的程序将产生该系统的阶跃响应模型, 并且绘制阶跃响应曲线（如图7.1所示）。 % 创建传递函数格式的模型 num = 1; den = 1 0.5; delt1 = 0.1; delay = 2; g = poly2tfd(num, den, delt1, delay);写成传递函数形式为 % 计算阶跃响应模型tfinal = 1.6;delt2 = delt1;nout = 1;plant = tfd2step(tfinal, delt2, nout, g); % 绘制阶跃响应曲线plotstep(plant)图 7.1 系统的阶跃响应曲线 我们也可以首先生成系

5、统的状态空间模型, 然后使用tf2ss函数和ss2step生成阶跃响应模型： num = 0 0 0 num; den = den 0 0; phi, gam, c, d = tf2ss(num, den); % 转换成状态空间形式 plant = ss2step(phi, gam, c, d, tfinal, delt1, delt2, nout); % 计算阶跃响应 通过MPC工具箱中的函数mpcinfo, 可以获取创建的阶跃响应模型中的信息： mpcinfo(plant) This is a matrix in MPC Step format. sampling time = 0.1 n

6、umber of inputs = 1 number of outputs = 1 number of step response coefficients = 16 All outputs are stable. 7.2.2 模型辨识 MATLAB中的MPC工具箱可以完成MISO系统模型的辨识。 基于输入v1(k), ， 和输出yl(k)的历史数据, 有(7.5)其中, 阶跃响应系数矩阵为 为了估计上述阶跃响应系数矩阵, 将SISO模型写成下面的形式其中, y(k)=y(k)-y(k-1), 脉冲响应系数hi=si-si-1。 如果输出端正在积分, 则(7.6) 其中, (y(k)=y(k)

7、-y(k-1), hi=hi-hi-1。 式（7.6）可以用来估计hi, 这样, hi和si可以写成 参数估计中通常将所有的变量缩放成相同的数量级。 这可以通过MPC函数autosc或scal完成。 整个系统数据可以写成下面的形式(7.7) 其中, Y包含了所有的系统输出信息（如果是稳定的, 则为y(k); 如果正在积分, 则为(y(k)）。 X包含了所有的输入信息（v(k)）。 包含所有被估计的参数（如果是稳定的， 则为hi; 如果正在积分, 则为hi）。 式（7.7）可以通过函数wrtreg得到。 参数可以利用多变量最小方差回归（mlr）或局部最小方差回归（plsr）计算得到。 下面的程序

8、（mpctutid）演示了具体的计算步骤： % 载入系统输入和输出数据, 这些数据由下面的传递函数和随机白噪声产生 % 从输入1到输出1的传递函数: g11=5.72exp(-14s)/ (60s+1) % 从输入2到输出1的传递函数: g12=1.52exp(-15s)/ (25s+1) load mlrdat ax, mx, stdx=autosc(x); mx=0*mx; % 将输入数据写成相同数量级的形式sx=scal(x, mx, stdx);n=35; % 一共具有35个系数xreg, yreg=wrtreg(sx, y, n); % 构造标准的辨识模型（见式（7.7）ninput

9、=2; plotopt=2;theta, yres=mlr(xreg, yreg, ninput, plotopt); % 计算脉冲响应系数theta=scal(theta, mx, stdx);% 将脉冲模型转换成阶跃模型, 用于MPC设计nout=1;delt=7;model=imp2step(delt, nout, theta);plotstep(model) % 绘制阶跃响应曲线图 7.2 实际输出与估计输出的误差分析曲线 7.2.3 无约束模型预测控制 MPC控制的基本原理可以用图7.4表示。 对于任何假定的现在和将来的控制输入u(k), u(k+1), ， u(k+m-1),系统将

10、来的输出y(k+1|k), y(k+2|k)， , y(k+p|k)都可以在一定的时间跨度p内进行预测。 从而可以计算m个现在和将来的控制输入(mp), 使得下面的二次指标最小 (7.8) 图 7.3 辨识后系统的阶跃响应曲线 图 7.4 MPC的控制过程 在缺省情况下, MPC将获得一个线性时不变的反馈控制规律 其中, Ep(k+1|k)是时间跨度为p的将来输出误差的预测向量, 它是在假定现在和将来的所有输入变量不发生变化（u(k)=u(k+1)=0）的情况下计算得到的。 对于开环稳定的系统, 闭环系统的稳定性是由KMPC决定的, 而KMPC与计算的时间跨度p、 m以及加权矩阵yl和ul有关

11、。我们一般无法知道关于这些参数保持系统闭环稳定的确切条件。 但是一个定性的描述是相对于p减小m可以增加系统的稳定度。 对于p=1, 闭环系统的稳定性可以由某个确定的m和时不变的输入输出权重加以保证。 通常, 权重矩阵ul用作可调常数, 因为增加ul一般可以使对系统的控制显得不那么“激进”, 从而增加系统的稳定性。 噪声滤波器的时间常数tfilter(1, :)和干扰时间常数tfilter(2, :)不会影响闭环系统的稳定性和对参考点变化或测量噪声的系统响应, 但是这些时间常数的设置会影响到系统的鲁棒性和对不可测干扰的响应。 增加噪声滤波器的时间常数可以使系统的鲁棒性增大, 使系统对不可测干扰的

12、敏感度降低。 而增加干扰的时间常数会使对系统的控制更加“激进”。 通过MPC工具箱所设计的控制器都可以无误差地跟踪阶跃信号（类型1）, 如果系统扰动模型或者系统自身具有积分特性, 则设计的控制器还可以无误差地跟踪斜坡信号。 下面的程序是MATLAB中MPC工具箱中演示例子（mpctutst.M）的一部分, 其中介绍了运用MPC工具箱进行模型预测控制系统设计的一般步骤和方法。 % 系统传递函数: g=5.72exp(-14s)/(60s+1)% 干扰的传递函数: gd=1.52exp(-15s)/ (25s+1)% 创建阶跃响应模型, 采样周期指定为7秒delt1=0;delay1=14;num

13、1=5.72;den1=60 1;g=poly2tfd(num1, den1, delt1, delay1);tfinal=245;delt2=7;nout1=1;plant=tfd2step(tfinal, delt2, nout1, g);delay2=15;num2=1.52;den2=25 1;gd=poly2tfd(num2, den2, delt1, delay2);delt2=7;nout2=1;dplant=tfd2step(tfinal, delt2, nout2, gd);% 计算MPC控制器的增益矩阵% 输出权重 = 1, 输入权重 = 0% 输入计算的时间跨度 = 5,

14、 输出计算的时间跨度 = 20model = plant; ywt=1; uwt=0;M=5; P=20;Kmpc1=mpccon(model, ywt, uwt, M, P);% 仿真系统dplant对不可测干扰和测量噪声的阶跃响应tend=245;r= ; usat= ; tfilter= ;dmodel= ;dstep=1;y1, u1=mpcsim(plant, model, Kmpc1, tend, r, usat, tfilter, . dplant, dmodel, dstep);dmodel=dplant; % 包含测量噪声y2, u2=mpcsim(plant, model,

15、 Kmpc1, tend, r, usat, tfilter, .dplant, dmodel, dstep);plotall(y1, y2, u1, u2, delt2); pause; % 对测量噪声干扰的阶跃响应, 暂停程序的运行, 观察计算的仿真曲线% 在另一组参数下重新计算MPC控制器的增益矩阵% 输出权重 = 1, 输入权重 = 10% 输入计算的时间跨度 = 5, 输出计算的时间跨度 = 20model=plant;ywt=1; uwt=10;M=5; P=20;mpc2=mpccon(model, ywt, uwt, M, P);% 仿真系统dplant对不可测干扰和测量噪声的

16、阶跃响应tend=245;r= ; usat= ; tfilter= ;dmodel= ;dstep=1;y3, u3=mpcsim(plant, model, Kmpc2, tend, r, usat, tfilter, .dplant, dmodel, dstep);dmodel=dplant; % 包含测量噪声y4, u4=mpcsim(plant, model, Kmpc2, tend, r, usat, tfilter, . dplant, dmodel, dstep);plotall(y3, y4, u3, u4, delt2);pause; % 暂停程序, 观察仿真结果 7.2.

17、4 闭环回路分析 除了直接进行仿真, MPC工具箱中还包含其它用来分析闭环系统稳定性和动态行为的函数和工具。 我们可以通过mpcc1函数获得闭环系统的状态空间描述, 然后可以利用smpcpole函数决定系统的极点位置。 下面的程序是mpctutst.M的一部分： % 构造无干扰的闭环系统, uwt = 0, 确定系统的极点 clmod = mpccl(plant, model, Kmpc1); poles = smpcpole(clmod); maxpole = max(poles) 计算结果为: maxpole = 1.0。 如果系统所有的极点均位于单位圆内或圆上, 则闭环系统是稳定的, 进

18、一步还可以分析闭环系统的频率响应。 对于多变量系统, 还可以利用svdfrsp函数计算系统奇异值关于频率的函数。 例如（mpctutst.M）% 计算灵敏函数和辅助灵敏函数的频率响应Cfreq = -3, 0, 30;ny = 1;in = 1:ny; % 辅助灵敏函数的输入为rout = 1:ny; % 辅助灵敏函数的输出为 ypfrsp, eyefrsp = mod2frsp(clmod, freq, out, in);plotfrsp(eyefrsp); % 计算灵敏度 pause;plotfrsp(frsp); % 计算辅助灵敏度pause; % 收有频率的幅值为1 7.2.5 有约束

19、的模型预测控制 MPC工具箱可以用来计算输入输出存在约束的系统。 假设系统输入变量约束(7.9) 输入速率约束 (7.10) 输出变量约束 (7.11) 如果系统存在上述的输入输出约束, 计算控制规律的每一步都需要求解二次方程, 最终所得的控制信号一般是非线性的, 这样系统的分析只能通过仿真进行。 让我们再看一下mpctutst.M程序中的相应处理（图7.5为计算后的仿真曲线）： % 分别在无约束和有约束情况下, 仿真系统dplant对系统阶跃干扰的响应 % 输出权重 = 1, 输入权重 = 0 % 输入计算的时间跨度 = 5, 输出计算的时间跨度 = 20图 7.5 有约束系统的仿真曲线 %

20、 输入信号的下限 = -0.4% 输入信号的上限 = inf% 输入的速率上限 = 0.1model = plant;ywt = 1; uwt = 0;M = 5; P = 20;tend = 245;r = 0;ulim = ;ylim = ; tfilter = ; dmodel = ;dstep = 1;y9, u9 = cmpc(plant, model, ywt, uwt, M, P, tend, r, .ulim, ylim, tfilter, dplant, dmodel, dstep);ulim = -0.4, inf, 0.1; % 加入约束y10, u10 = cmpc(p

21、lant, model, ywt, uwt, M, P, tend, r, .ulim, ylim, tfilter, dplant, dmodel, dstep);plotall(y9, y10, u9, u10, delt2);7.3 基于状态空间模型的模型预测控制 7.3.1 MPC表示的状态空间模型 考虑图7.6所示的系统对象框图。 MPC工具箱中使用的离散线性时不变（LTI）系统状态方程的一般描述为(7.12)图 7.6 MPC被控系统的框图 其中, x表示n个状态变量的状态向量, u表示nu个输入变量, d代表nd个可测的但只有变化的输入（也就是可测的干扰信号）, w代表nw个不可

22、测干扰。 y是ny个系统输出, z是相应的测量噪声, 和u等是一定大小的常数矩阵, 变量y(k)表示没有加入测量噪声的系统输出。 定义 在MPC工具箱中, 系统的状态方程模型是以特殊的格式保存的, 称为mod格式。它是一个包含系统状态空间基本矩阵、 、 C、 D和其它一些信息的矩阵。 MPC工具箱提供了一系列函数完成其它模型向mod格式的转换。 1） SISO连续时间传递函数的转换 假设连续时间传递函数的一般表示为(7.13) 例 7.1 考虑下面的传递函数模型, 创建该模型的MPC mod格式表示。 解： 具体的代码如下： G = poly2tfd(-13.6 1, 54.3 113.5 1

23、, 0, 5.3); % 创建传递函数模型mod = tfd2mod(2.1, 1, G); % 将传递函数模型转换成MPC mod格式, 采样周期为2.1秒 2） SISO离散时间传递函数的转换 假设离散时间传递函数的一般表示为(7.14) 其中的Td表示系统的时间延迟。 使用tf或poly2tfd函数可以创建上述模型的系统。 然后通过tfd2mod命令可以将传递函数形式转换成MPC表示的mod格式。 例 7.2 将下面的离散SISO系统转换成MPC mod格式 解： 具体代码为 G = poly2tfd(-0.1048 0.1215 0.0033, 1 -0.9882 0.0082, 2.

24、1, 3);mod = tfd2mod(2.1, 1, G); 注意， 这里的poly2tfd和tfd2mod命令指定了相同的采样周期（delta = 2.1）。 实际上, 也可以在两个命令中指定不同的采样周期。 3） MIMO传递函数的转换 假定MIMO系统具有下面的传递函数矩阵描述 其中, gi,j表示第j个输入到第i个输出之间的传递函数。 如果所有ny个输出都可以测量, 所有nu个输入都是控制变量, 则直接调用tfd2mod命令就可以得到正确的mod格式模型。 例如, 对于下面的三输入二输出系统(7.15)下面的命令将该系统转换成mod格式(其中采样周期为T=4 s)： g11 = po

25、ly2tfd(12.8, 16.7 1, 0, 1); g21 = poly2tfd(6.6, 10.9 1, 0, 7); g12 = poly2tfd(-18.9, 21.0 1, 0, 3); g22 = poly2tfd(-19.4, 14.4 1, 0, 3); g13 = poly2tfd(3.8, 14.9 1, 0, 8); g23 = poly2tfd(4.9, 13.2 1, 0, 3); pmod = tfd2mod(4, 2, g11, g21, g12, g22, g13, g23); 如果假设其中的第三个输入是不可测的干扰, 相应的系统描述变成(7.16) 在这种情

26、况下, 需要通过分别指定控制变量、 可测干扰和不可测干扰的输入数目来重载tfd2mod命令的缺省模式。 通常只要修改所得模型的第1行的相应内容即可。 其中的具体含义如下： 第1列元素T： 系统的采样周期。第2列元素n： 系统的状态个数。第3列元素nu： 控制输入变量的数目。第4列元素nd： 测量干扰的数目。第5列元素nw： 不可测干扰的数目。第6列元素nym： 可测量输出的数目。第7列元素nyu： 不可测量输出的数目。 4） 离散或离散状态方程的转换 使用c2dmp（将连续系统转换成离散状态方程形式）和ss2mod（将离散状态方程写成mod格式）命令可以将连续时间状态方程模型转换成MPC mo

27、d格式的模型。 当然, 如果已经得到离散状态方程形式, 第一步就不需要了。 例如, 假定a、 b、 c和d是描述一个连续系统的矩阵, 为了将它转换成采样周期为T=1.5 s的mod格式, 可以输入下面的命令: phi, gam = c2dmp(a, b, 1.5); mod = ss2mod(phi, gam, c, d, 1.5); 5） mod格式模型的组合 MPC工具箱中的addmod、 addumd、 appmod、 paramod和sermod等命令允许用户将多个简单的mod格式模型组合成更为复杂的系统结构。 例如, addmod可以将多个可测扰动的影响集中到单个系统中, 如图7.7

28、所示。 图 7.7 多个mod模型的组合 而pmod可以集中多个控制变量u(k)的影响, 也可以对输出变量y(k)上的不可测干扰w(k)进行类似的处理。 dmod命令可以组合同一输出变量上的多个测量干扰。 一旦完成这些信号组合的工作后, 就可以使用addmd命令产生组合后的最终模型： model = addmd(pmod, dmod); 6） mod格式向其它类型模型的转换 函数mod2ss可以将mod格式的模型转换成标准的离散状态方程形式： phi, gam, c, d, minfo = mod2ss(mod); 其中, phi、 gam、 c和d分别是下面状态方程的系数矩阵。 x(k+1)

29、=x(k)+u(k) y(k)=Cx(k)+Du(k)向量minfo则包含mod模型中的第一行中的内容。 7.3.2 基于状态空间模型的无约束MPC设计 一旦得到被研究系统的mod格式描述, 就可以使用下面的命令完成系统设计的工作： smpccon： 计算无约束控制器的增益矩阵。 smpcest： 设计状态观测器（可选）。 smpcsim： 计算闭环系统在指定输入信号下的响应。 plotall(ploteach)： 绘制计算得到的仿真曲线。 另外, 还可以使用下面的命令分析闭环系统的特定属性： smpcc1： 产生闭环系统模型（系统对象加上控制器)。 smpcgain： 计算闭环系统的增益矩阵

30、。smpcpole： 计算闭环系统的极点。mod2frsp(plotfrsp)： 计算并绘制闭环系统频率响应曲线。 svd2frsp： 计算频率响应的奇异值。 例 7.3 对于式（7.16）描述的MIMO系统, 来设计该系统的MPC控制器。 解： 首先将系统定义成mod格式, 下面的程序将完成这一工作（采样周期设置为2秒） delt = 2; ny = 2; g11 = poly2tfd(12.8, 16.7 1, 0, 1); g21 = poly2tfd(6.6, 10.9 1, 0, 7); g12 = poly2tfd(-18.9, 21.0 1, 0, 3); g22 = poly2

31、tfd(-19.4, 14.4 1, 0, 3);umod = tfd2mod(delt, ny, g11, g21, g12, g22);% 定义u 输入的效果g13 = poly2tfd(3.8, 14.9 1, 0, 8);g23 = poly2tfd(4.9, 13.2 1, 0, 3);dmod = tfd2mod(delt, ny, g13, g23);% 定义w 输入的效果pmod = addumd(umod, dmod); % 将上面得到的两个模型进行组合 接下来设计系统的无约束MPC控制器。 设计参数本质上应与基于阶跃响应模型的函数（参考前一节的内容）相同。 在本例中, 选择

32、如下的设计参数: imod = pmod; % 假设建立的模型是准确的 ywt = ; % 两个输出采用缺省（单位）的权重 uwt = ; % 两个输入采用缺省（零）的权重 P = 5; % 预测时间宽度 M = P; % 控制的时间宽度 Ks = smpccon(imod, ywt, uwt, M, P); 下面通过仿真输出y1的阶跃响应检查所设计的控制器。 图7.8是仿真的结果。 tend=30; % 仿真的时间长度 r = 1 0; % 两个输出的设置点 y, u = smpcsim(pmod, imod, Ks, tend, r); plotall(y, u, delt)图 7.8 输出阶跃响应的仿真曲线 7.3.3 基于状态空间模型的有约束MPC设计 函数scmpc可以处理控制变量和输出变量上的不等式约束问题。 解决有约束状态空间MPC问题的基本步骤是首先使用前一小节讨论的函数, 使用系统的状态空间模型分析无约束情况下的预测时间跨度、 控制时间跨度、 系统的输入输出权重和系统的状态观测器等等。 然后再考虑有约束的情况, 并利用scmpc函数求解该问题。 下面的例子将演示求解有约束MPC问题的基本过程。 让我们再看