导航

Global optimizationnonconvex quadratic programmingnonconvex polynomial programminggeneral nonconvex programmingnonconvex semidefinite programmingReferences

Global optimization

内置求解器BMIBNB可以求出全局最优解，并且在规模较小的问题上表现出很好的鲁棒性，但是求解速度较慢.

nonconvex quadratic programming

第一个例子具有一个concave quadratic constraint，在分支定界中，出现三种不同优化问题：Upper bounds 使用bmibnb.uppersolver，Lower bounds使用bmibnb.lowersolver，bound tightening使用bmibnb.lpsolver.

clear all; x1 = sdpvar(1, 1); x2 = sdpvar(1, 1); x3 = sdpvar(1, 1); p = -2*x1+x2-x3; % 设置约束 F = [x1*(4*x1-4*x2+4*x3-20)+x2*(2*x2-2*x3+9)+x3*(2*x3-13)+24>=0, 4-(x1+x2+x3)>=0, 6-(3*x2+x3)>=0, x1>=0, 2-x1>=0, x2>=0, x3>=0, 3-x3>=0]; ops = sdpsettings('solver', 'bmibnb'); optimize(F, p, ops)

第二个例子是一个nonconvex quadratic programming

clear all x1 = sdpvar(1,1); x2 = sdpvar(1,1); x3 = sdpvar(1,1); x4 = sdpvar(1,1); x5 = sdpvar(1,1); x6 = sdpvar(1,1); p = -25*(x1-2)^2-(x2-2)^2-(x3-1)^2-(x4-4)^2-(x5-1)^2-(x6-4)^2; F = [(x3-3)^2+x4>=4,(x5-3)^2+x6>=4,x1-3*x2<=2, -x1+x2<=2, x1-3*x2<=2, x1+x2>=2,6>=x1+x2>=2,1<=x3<=5, 0<=x4<=6, 1<=x5<=5, 0<=x6<=10, x1>=0,x2>=0]; options = sdpsettings('solver','bmibnb'); optimize(F,p,options)

发现当前目标值为-313，而下界值为-310，需要缩小目标值之间的gap.

The hard part is often not finding good solutions, but proving their optimality.

options=sdpsettings('solver', 'bmibnb', 'bmibnb.relgaptol', 1e-4); optimize(F, p, options)

nonconvex polynomial programming

多项式规划(polynomial programs)一般会转为双线性规划(bilinear programs)

sdpvar x y F = [x^3+y^5<=5, y>=0]; ops = sdpsettings('verbose', 1, 'solver', 'bmibnb'); optimize(F, -x, ops)

general nonconvex programming

sdpvar x y p = sin(1+y*x)^2+cos(y*x); optimize([-1 <= [x y] <= 1],p,sdpsettings('solver','bmibnb'));

nonconvex semidefinite programming

The following problem is BMI, bilinear matrix inequality problem.

yalmip('clear') x = sdpvar(1,1); y = sdpvar(1,1); t = sdpvar(1,1); A0 = [-10 -0.5 -2;-0.5 4.5 0;-2 0 0]; A1 = [9 0.5 0;0.5 0 -3 ; 0 -3 -1]; A2 = [-1.8 -0.1 -0.4;-0.1 1.2 -1;-0.4 -1 0]; K12 = [0 0 2;0 -5.5 3;2 3 0]; F = [x>=-0.5, x<=2, y>=-3, y<=7]; F = [F, A0+x*A1+y*A2+x*y*K12-t*eye(3)<=0]; options = sdpsettings('solver','bmibnb'); optimize(F,t,options);

Starting from release R20200930, the default behavior to attack BMIs in BMIBNB is by employing a cutting plane strategy for the upper bound generation.

第二个BMI问题要求求解一个Linear Qudratic Control问题(LQ，线性二次型问题)，需要找到约束条件下的状态反馈矩阵(find a state-feedback matrix with bounded elements).

For the global code to work, global lower and upper and bound on all complicating variables (involved in nonlinear terms) must be supplied, either explicitly or implicitly in the linear constraints.

本例中，变量$\mathbf{K}$在原问题中是受约束的条件，但是$\mathbf{P}$需要手动设置约束

%% LQR yalmip('clear') A = [-1 2; -3 -4]; B = [1; 1]; P = sdpvar(2, 2); K = sdpvar(1, 2); F = [P>=0, (A+B*K)'*P+P*(A+B*K) <= -eye(2)-K'*K]; F = [F, -0.1<=K<=0.1]; ops = sdpsettings('solver', 'bmibnb'); optimize(F, trace(P), ops)

注意到BMIBNB给出了无界变量的警告，容易推测出$\mathbf{P}$的对角元素为非负，但是没有对元素上界进行限制，所以求解器无法推测出关于非对角元素的限制，由于缺少限制条件，BMIBNB可能无法收敛.

Notice that we have some degree-of-freedom in how we model this problem. The Lyapunov function which is nonlinear in $K$ and $P$ can be paritally linearized by applying a Schur Complement.

F = [P>=0, [-eye(2) - ((A+B*K)'*P+P*(A+B*K)) K';K 1] >= 0]; F = [F, K<=0.1, K>=-0.1]; %加入限制条件 ops = sdpsettings('solver','bmibnb'); optimize(F,trace(P),ops)

求解一个衰减速率最大化问题，该BMI问题的求解效率就会变得很低

%% yalmip('clear'); A = [-1 2;-3 -4]; t = sdpvar(1,1); P = sdpvar(2,2); F = [P>=eye(2), A'*P+P*A <= -2*t*P]; F = [F, t >= 0]; ops = sdpsettings('solver','bmibnb'); optimize(F,-t,ops)

For this particular problem, the reason is easy to find. The original BMI is homogeneous, and to guarantee a somewhat reasonable solution, we artificially added the constraint $P⪰I$ instead of $P\succ 0$.

%% F = [P>=0, trace(P)==1, A'*P+P*A <= -2*t*P]; F = [F, t >= 0]; optimize(F,-t,ops)

For this problem, we can easily find more interesting cutting planes. The decay-rate BMI together with the constant trace implies $trace(A^{T}P+PA)\le -2t$. Adding this redundant cut leads to a finite lower bound.

%% F = [P>=0, trace(P)==1, A'*P+P*A <= -2*t*P]; F = [F, t >= 0]; F = [F, trace(A'*P+P*A)<=-2*t]; % add cut plane optimize(F,-t,ops);

A Schur complement on the decay-rate BMI gives us yet another linear SDP cut which improves the node relaxation even more.

%% F = [P>=0,A'*P+P*A <= -2*t*P, t >= 0]; F = [F, trace(P)==1]; F = [F, trace(A'*P+P*A)<=-2*t]; F = [F, [-A'*P-P*A P*t;P*t P*t/2] >= 0]; % Schur complement optimize(F,-t,ops);

有效的cut与效率问题

By adding valid cuts, the relaxations are possibly tighter, leading to better bounds. A problem however is that we add additional burden to the local nonlinear solver used for the upper bounds.The additional cuts are redundant for the local solver, and most likely detoriate the performance. To avoid this, cuts can be explicitly specified by using the command cut.

F = [P>=0,A'*P+P*A <= -2*t*P,100 >= t >= 0]; F = [F, trace(P)==1]; F = [F, trace(A'*P+P*A)<=-2*t]; F = [F, cut([-A'*P-P*A P*t;P*t P*t/2]>=0)]; % cut指令显式指定 optimize(F,-t,options);

References

Global optimization