Navigator

mutiparametric programminggeneric exampleSimple MPC exampleMixed integer multiparametric programming

mutiparametric programming

多参数规划求解需要安装MPT求解器 MPT安装脚本下载

generic example

以下的二次规划(QP)，决策变量为$\mathbf{z}$，参数向量为$\mathbf{x}$

%% A = randn(15,3); b = rand(15,1); E = randn(15,2); z = sdpvar(3,1); x = [0.1;0.2]; F = [A*z <= b+E*x]; obj = (z-1)'*(z-1); sol = optimize(F,obj); value(z)

为了得到关于$\mathbf{x}$的参数解，需要调用函数solvemp，并且显式告知求解器$\mathbf{x}$是参数变量，并且为参数变量定义范围

%% x = sdpvar(2, 1); F = [A*z<=b+E*x, -1<=x<=1]; sol = solvemp(F, obj, [], x)

输出值是MPT的解结构，设置参数得到对应参数下的解

% 设置参数，得到对应解 xx = [0.1; 0.2]; [i, j]=isinside(sol{1}.Pn, xx); sol{1}.Fi{j}*xx+sol{1}.Gi{j}

可以利用YALMIP提供的nonlinear operator framework整合流程

%% [sol, diagnostics, aux, value_func, optimal_z]=solvemp(F, obj, [], x); % assign and value assign(x, [0.1; 0.2]); value(optimal_z)

MPT重载了绘图函数，可以绘制出不同参数值下的最优解情况

% plot plot(value_func); plot(optimal_z);

Simple MPC example

定义线性系统的数值，预测矩阵和当前状态的变量$x$，控制序列为$U(x)$，对于一个5阶段的MPC系统描述如下（create_CHS是一个cheat_function用于产生当前状态$x$和输入控制序列$U$的数值矩阵）

N = 5; A = [2 -1; 1 0]; B = [1; 0]; C = [0.5 0.5]; [H, S]=create_CHS(A, B, C, N); x = sdpvar(2, 1); U = sdpvar(N, 1);

预测输出为当前状态和控制序列的函数

Y = H*x+S*U;

设置费用函数为输出序列和控制序列消耗的二次项之和

obj = Y'*Y+U'*U;

设置约束条件，末端约束和控制向量约束

F = [1>=U>=-1, 1>=Y(N)>=1];

在$|x|\le 5$条件下求解控制向量$U(x)$

The explicit solution $U(x)$ is obtained by calling solvemp with the parametric variable $x$ as the fourth argument. Additionally, since we only are interested in the first element of the solution $U(x)$, we use a fifth input to communicated this.

F = [F, 5>=x>=-5]; % solve [sol, digno, aux, value_func, Optimizer]=solvemp(F, obj, [], x, U(1));

求出控制序列如图所示 价值函数如图

Mixed integer multiparametric programming

通过添加二元变量(binary variable)可以求解扩展MPC问题

N = 3; A = [2 -1;1 0]; B = [1;0]; C = [0.5 0.5]; [H,S] = create_CHS(A,B,C,N); x = sdpvar(2,1); U = sdpvar(N,1); Y = H*x+S*U; objective = norm(Y,1) + norm(U,1); F = [1 >= U >= -1]; F = [F, 5 >= x >= -5];

加入额外的逻辑约束后，这类非凸算子引入许多二元变量，此时MPC问题比较有效率的求解方式是动态规划

F = [F, ismember(U, [-1:1/3:1])]; [sol,diagnostics,aux,Valuefunction,Optimizer] = solvemp(F,objective,[],x,U(1)); plot(Optimizer);