MATLAB约束最优化之罚函数法、障碍函数法和SQP方法

1. 罚函数法

罚函数方法包括外点法和内点法。外点法又叫外罚函数法，顾名思义，迭代点再约束条件的可行域之外，既用于不等式约束又可用于等式约束。

同样地，罚函数方法又叫序列无约束最小化方法。顾名思义，罚函数方法是通过构造惩罚项来消除约束条件，从而转化为无约束最优化问题。

Matlab通过fminunc来求解无约束最优化问题。下述程序把整个程序封装到函数runfminunc里面，返回最优解xsol和最优值fval、迭代过程值history和搜索路径searchdir。

其中函数fminunc是求解无约束最优化问题的，参数@objfun是调用目标函数句柄，x0是初始值，options是优化参数，其中@fminunc是优化方法，@outfun定义了迭代过程中每一步怎么做，这里是每一步都保存当前值到结构体history中，并保存迭代中的搜索方向searchdir，Display设为iter即可，Algorithm表示用quasi-newton拟牛顿法求解。

[xsol,fval,history,searchdir] = runfminunc

function [xsol,fval,history,searchdir] = runfminunc
history.x = [];
history.fval = [];
searchdir = [];
x0 = [0, 0];
options = optimoptions(@fminunc,'OutputFcn',@outfun,...
    'Display','iter','Algorithm','quasi-newton');
[xsol,fval] = fminunc(@objfun,x0,options);

    function stop = outfun(x,optimValues,state)
        stop = false;
        
        switch state
            case 'init'
                disp('start');
            case 'iter'
                history.fval = [history.fval; optimValues.fval];
                history.x = [history.x; x];
                searchdir = [searchdir;...
                    optimValues.searchdirection'];
            case 'done'
                disp('end');
            otherwise
        end
    end

    function f = objfun(x)
        M = 5000;
        f = (x(1)-2.25)^2+(x(2)-2)^2;
        g = f+M*max(x(1)^2-x(2),0)+M*max(x(1)+x(2)-6,0)-M*min(x(1),0)-M*min(x(2),0);
    end

end

2. 障碍函数法（又叫内点法、内罚函数法）

障碍函数法又叫内点法，也是罚函数法的范畴，内点法一般只适用于 不等式约束 的优化问题。

障碍函数：

Example：

min $f(x_1,x_2)=-x_1x_2$, s.t. $1-x_2^2-x_2^2 \geq 0$。

对该问题构造障碍函数为：$B(x)=-\frac{1}{\beta} log(1-x_1^2-x_2^2)$

障碍因子$\beta=1 \rightarrow \infty$。

内点法就要用fmincon方法来求解了，在MATLAB中，fmincon方法如下：

其中Algorithm字段设置为了‘interior-point’，表示通过内点法求解fmincon所定义的约束优化问题，其实fmincon的options中缺省了Algorithm字段的话，默认值就是内点法。除此之外还有‘sqp’等。

options = optimoptions(@fmincon,'OutputFcn',@outfun,... 
    'Display','iter','Algorithm','interior-point');
[xsol,fval] = fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options);

完整代码如下：

其中fmincon方法的参数按顺序分别是：

目标函数、初始解、线性等式约束A和b、线性不等式约束Aeq和beq、解的范围lb和ub，这里一律设置为了空‘[]’，表示没有。

@confun是定义了非线性约束的句柄。

function [xsol,fval,history,searchdir] = interiorfmincon
 
% Set up shared variables with outfun
history.x = [];
history.fval = [];
searchdir = [];
 
% Call optimization
x0 = [0, 0];
options = optimoptions(@fmincon,'OutputFcn',@outfun,... 
    'Display','iter','Algorithm','interior-point');
[xsol,fval] = fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options);

 function stop = outfun(x,optimValues,state)
     stop = false;
 
     switch state
         case 'init'
             disp('start');
         case 'iter'
           history.fval = [history.fval; optimValues.fval];
           history.x = [history.x; x];
         case 'done'
             disp('end');
         otherwise
     end
 end
 
 function f = objfun(x)
     f = (x(1)-2.25)*(x(1)-2.25) + (x(2)-2)*(x(2)-2);
 end
 
 function [c, ceq] = confun(x)
     % Nonlinear inequality constraints
     c = [x(1)*x(1)-x(2);
         x(1)+x(2) - 6];
     % Nonlinear equality constraints
     ceq = [];
 end
end

调用[xsol, fval, history, searchdir] = iteriorrunfmincon即可

3. SQP方法（即序列二次规划法）

SQP方法把约束最优化问题转换为二次规划，迭代求解。代码和内点法几乎一模一样，只有一个地方不同，如下：

options = optimoptions(@fmincon,'OutputFcn',@outfun,... 
    'Display','iter','Algorithm','sqp');
[xsol,fval] = fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options);