基于SMPC的H∞优化问题的求解的方法

1、标题：基于SMPC 的H优化问题的求解的方法时间：2016-4-17内容：1.首先介绍一个简单的例子的MLD系统的推导方法2.分析SMPC框架下的优化问题的线性化方法3.根据系统模型以及线性化问题定义，在SMPC框架下推导得到具有场景信息的LP求解参数目录1.例子：32问题：43.模型线性化、约束线性化61.例子： 那么，关于条件 ，引入二进制变量 则有 结合 ，令 ，那么 通过逻辑变换规则，那么（2）式可以等效为 那么（1）式可以写成 设 ， 那么(3-4)式可以进一步等效为 那么(5)式可以进一步等效为 （6-7）可以被概括为MLD系统 其中， 其中，连续变量的维数为 ， 集合维数为 2问

2、题：在预测控制中，普遍采用二次型的优化性能指标，它对应于预测偏差的二范数，反应了对于输出值良好跟踪期望值的性能要求。当控制要求侧重于良好的镇定性能，即所有的输出值都保持在设定值附近，不希望出现过大的偏差。用以刻画控制要求的更直接的数学度量是要极小化最有可能出现的最大偏差（min-max）,即 相当于将二次型性能指标推广到无情范数形式 其中，向量 的无穷范数定义为 当（11）中的无穷范数对空间（全部输出量）和时间（整个优化时域）双重取值时，其具体表达式就是（10），注意到，性能指标（10）具有线性形式，如果对应的模型和物理约束也是线性的，那么这一优化问题可以用线性规划方式求解。即，当满足一下三点

3、时可以通过LP进行求解：一、 性能指标线性化二、 模型线性化三、 约束线性化一、 性能指标线性化 （10）式涉及到极值和结对子项，与标准LP问题形式尚存在差距，需要将其转化成标准LP问题。根据Morari M. Robust Model predictive Control.介绍的方法：对于m个输入，p个输出的对象模型，预测时域长度为N。设从k时刻起 个控制量 作用下对未来 步的输出闭环预测为 其中， 那么（13）式可以等价为 设k时刻的滚动优化问题是对输入输出有约束情况下的无穷范数有福问题： 引入 那么（16）可等效为 其中， 是 的第 个分量 为了将（10）式转变成线性规划问题，定义 显然

4、，这一最优 值是满足 那么（18）可以等效为 其中， 是 的线性形式。LP的标准形式为： 3.模型线性化、约束线性化考虑离散线性离散时间系统 根据例子的方法，引入0-1变量 可将系统转变成，MLD形式的 其中，优化问题的无穷范数可以定义为 根据在例子中介绍的代换方法： 可以对（24）的约束可以作同样的变换：实现代码:MLD.A = S.A;MLD.B1 = S.B;MLD.Bw = S.D*S.eps;MLD.C = S.C;其中，S.A、S.B、S.C、S.D：是系统模型MLD.E1 = zeros(size(S.Ax,1),nu);-S.Au;MLD.E4 = S.Ax;zeros(siz

5、e(S.Au,1),nx);MLD.E5 = S.bx;S.bu;MLD.Ew = zeros(size(MLD.E1,1),nw); ， ， 注意：程序中的 与（24）式约束相符合，即 其中，（28）式与文章“Control of systems integrating logic, dynamics,and constraints1 Alberto Bemporad, Manfred Morari*”的（10）式相符合；注意到：在程序中，并没有考虑 ， 定义为空。接下来获取LP求解过程参数：实现代码:nu = size(B1,2);nd = size(B2,2); %空nz = size(

6、B3,2); %空nx = size(A1,1);ne = size(E1,1);% Slack variablesnq = N+N; % Here we're assuming nu>0, nx>0, nz>0 % u,delta variablesno = (nu)*N; % total number of varsnvar = nq+no; % optim vector: eps_u,eps_x,eps_z,u,d,z % Ô¼ÊøÌõ¼þÊýnc = N*(ne

7、+2*nu+2*nx);根据线性规划问题定义： 求解 那么需要得到，基于场景的矩阵参数为 实现代码：for i = 1:rC = smpc_bsmpc_sce_lp(sys,L,treeinf.Si);%构建每条场景的结构C. f = f;treeinf.Pi*C.f; b = b;C.b; Cx = Cx;C.Cx; Cru = Cru;C.Cr.u; Crx = Crx;C.Cr.x; A(i-1)*nq + 1):(i-1)*nq + nq),(i-1)*nz + 1):(i-1)*nz + nz) = C.A;end下面将分别介绍以上7个主要的参数矩阵的构建方式：l 常参数部分：1.

8、C.f实现代码：% Column array containing the objective function coefficients.smpc.f=ones(1,nq) zeros(1,no)'注意到： 其中， 之所以为 维，是因为我们的问题的输入 是二维的。% Slack variablesnq = N+N; %假设 nu>0, nx>0, nz>0.% u,delta variablesno = (nu)*N;2. C.Cr.u % Epsilon-u inequalitiesfor i=1:N, （注意for循环的目的是将树的动态信息体现到矩阵中去） Cr

9、=Cr;-Qu; Cr=Cr;Qu;endsmpc.Cr.u(1:N*2*nu,:)=Cr;由于 C.Cr.u 是输入参考值的系数矩阵，而在目标函数中，输入的阶段成本函数为 ，因此，去绝对值后 的系数就可以确定了，且仅与 Qu 有关为什么存在正负两部分组成，见（31）式3. C.Cr.x % Epsilon-x inequalitiesfor i=1:N, Cr=Cr;-Qx; Cr=Cr;Qx;endsmpc.Cr.x(2*N*nu+1:2*N*nu+2*N*nx,:)=Cr;同2l 不等式参数部分：注意到：不等式参数部分有三部分组成：1) Epsilon-u inequalities：对应

10、于输入物理约束不等式2) Epsilon-x inequalities：对应于状态物理约束不等式3) MLD inequalities：对应于MLD系统约束不等式4. C.bsmpc.b=Sp;%不依赖于x(t),r(t)的约束系数矩阵。【常量】实现代码：（注意for循环的目的是将树的动态信息体现到矩阵中去）for i=1:N, Sp=Sp;zeros(nu,1); % Epsilon-u inequalitiesSp=Sp;zeros(nu,1); Sp=Sp;zeros(nx,1); % Epsilon-x inequalities Sp=Sp;zeros(nx,1); Sp=Sp; E5

11、; % MLD inequalitiesend为什么会存在三部分呢？在我们考虑的问题中，我们既有状态约束 ，也有输出约束 。参照例子中讲到的变换方式，它们分别引入 和 （0-1）变量，二者变换得到的不等式就分别称作：Epsilon-x inequalities 和 Epsilon-u inequalities 。根据，描述的 MLD inequalitie是华为MLD形式后的不等书与 和 的关系5. C.Asmpc.A=Gp;for i=1:N,% Epsilon-u inequalitiesGp=Gp;zeros(nu,i-1) -ones(nu,1) zeros(nu,(nq-i+nu*(

12、i-1) -Qu zeros(nu,no-nu*i);注：Gp=Gp;zeros(nu,i-1) -ones(nu,1) zeros(nu,(nq-i+nu*(i-1) Qu zeros(nu,no-nu*i);注：其中 注意到，为什么在 不等式的变量系数矩阵中会存在这样正负的关系结构。是因为 不等式表示的是物理约束 % Epsilon-x inequalities Gpxu=; Gpxd=; Gpxz=; Gpxw=;for j=1:i aux = eye(size(A1); for k = 1:(i-j) aux = aux*Ai-k+1; end Gpxu=Gpxu -Qx*aux*Bj

13、; Gpxw=Gpxw -Qx*aux*Bw; end Gpx=Gpxu zeros(nx,(nu*N)-size(Gpxu,2) Gpxd zeros(nx,(nd*N)-size(Gpxd,2) Gpxz zeros(nx,(nz*N)-size(Gpxz,2); Gp=Gp;zeros(nx,N+i-1) -ones(nx,1) zeros(nx,N-i) Gpx;aux = eye(size(A1); for k = 1:i aux = Ak*aux; endGp=Gp;zeros(nx,N+i-1) -ones(nx,1) zeros(nx,N-i) -Gpx;推导过程同上：Gpeu

14、=; Gped=; Gpez=; Gpew=;% MLD inequalities for j=1:i-1 aux = eye(size(A1); for k = 1:(i-j-1) aux = Ai-k*aux; end aux = -E4*aux; Gpeu=Gpeu aux*Bj; Gpew=Gpew aux*Bw; end Gpeu=Gpeu -E1; Gpew=-Gpew Ew; Gpe=Gpeu zeros(ne,nu*N-size(Gpeu,2) Gped zeros(ne,nd*N-size(Gped,2) Gpez zeros(ne,nz*N-size(Gpez,2);Gp=Gp; zeros(ne,nq) Gpe;推导过程同上：end6. C.Cx（同5）smpc.Cx=Fp;for i=1:N, Fp=Fp;zeros(nu,nx); % Epsilon-u inequalitiesFp=Fp;zeros(nu,nx); aux = eye(siz