active-set 的SQP算法

%-------------------------------------------------------------------------- % % 子函数 GQP % % 采用有效集算法求解一般严格凸二次规划: % // Generalized strictly convex Quadratic Programming % min f(x)=0.5*x'Gx+C'x % s.t. A*x>=b % % Author : ZDL Date : 5/13/2011 % %-------------------------------------------------------------------------- % % Input Variables : % G - n阶正定对称矩阵 % C - n*1列向量 % A - 列满秩矩阵 % b - 列向量 % epsilon - 迭代控制误差 % % Output Variables : % x - 最优解，列向量 % lambda - Lagrange乘子列向量 % % Local Variables : % cs - 约束对应的结构数组 % d - 搜索方向 % alpha - 步长 % % Invoking Functions : % init_feasible - 产生初始可行点函数 % AssignStructValue - 填充约束结构数组函数 % detach - 结构数组域分离函数 % QP - 求解等式约束QP函数 % attach - 有效约束对应的Lagrange乘子与 % 约束结构的lambda字段相关联函数 % %-------------------------------------------------------------------------- function [x,lambda]=GQP(G,C,A,b,epsilon) % Find the minimum of the generalized strictly convex quadratic % programming using active set algorithm. % 计算初始可行点 //线性规划问题可能无解 try x0=init_feasible(A,b); catch exception [m,n]=size(A); x0=zeros(n,1); end % 约束结构数组cs 初始化 cs=AssignStructValue(x0,A,b); % 迭代点初始化 xk=x0; % 采用有效集算法迭代求解QP问题 num=0; while 1 num=num+1; % 分离出有效约束/紧约束 [AA1,bb1,index1]=detach(cs,true); % 分离出松约束 [AA2,bb2,index2]=detach(cs,false); % 构造只具有等式约束的QP问题 gk=G*xk+C; bb=zeros(size(bb1)); % 求解只具有等式约束的QP问题，给出当前步下降可行搜索方向d 以及相应的 % Lagrange乘子向量lamb [d,lamb]=QP(G,gk,AA1,bb); % 判断 d=0，作为迭代停止准则 if norm(d)<epsilon || num>50 if ~isempty(lamb) % 计算不等式约束集合中有效约束对应的最小Lagrange乘子以及相应的位置 [L,s]=min(lamb); if L<0 % 若 d=0 且Lagrange乘子向量有负分量，则从有效约束集中删除第s个 % 不等式约束（只是在当前点为有效约束），重新求解QP问题 cs(index1(s)).ActiveSet=false; continue; end end % 若 d=0 且不等式约束对应的Lagrange乘子向量非负，或为空集，则输出最优解 x=xk; % 计算最优值点对应的约束结构cs cs=AssignStructValue(x,A,b); % 关联有效约束对应的Lagrange乘子 cs=attach(cs,lamb,index1); %输出原问题约束对应的Lagrange乘子向量 lambda=[cs.lambda]'; %迭代结束并返回 break; else % 若 d>0，则沿d方向作线性搜索，寻求更好点 Alpha_k=[]; % 存储到达各松约束边界的最大有限步长 T=[]; % 存储具有最大有限步长的松约束对应的编号 NUM=length(index2); % 计算松约束个数 I=1; % 在所有松约束集合中查找具有最大有限步长的松约束并记录相应的步长 for i=1:NUM % 只有 AA2(i,:)*d<0 时对步长有限制 if AA2(i,:)*d<0 Alpha_k(I)=(bb2(i)-AA2(i,:)*xk)/(AA2(i,:)*d); T(I)=index2(i); I=I+1; end end t=0; % 记录当前搜索方向最先到达的松约束编号 alpha_k=2; % 记录当前搜索方向最先到达的松约束对应的步长 if NUM>0 % 计算当前搜索方向最先到达的松约束步长以及编号 [alpha_k,tt]=min(Alpha_k); t=T(tt); end % 判断当前搜索方向下最优值点是否落在边界上 alpha=min([1,alpha_k]); % 求解步长 xk=xk+alpha*d; % 更新迭代点 % 若当前搜索方向下最优值点落在边界上，则相应的松约束成为有效约束，原来的 % 有效约束仍保持不变 if alpha<1 cs(t).ActiveSet=true; % 添加有效约束 end end end