非线性规划MATLAB代码

    下面是三个非线性规划领域的算法。课堂上给予了详细的讲解，在实践环节让学生编程实现，从而可以实验复杂一些的例子，加深对算法的理解。下面共有四个程序grad，simplelinesearch，bfgs和phr，全部使用MATLAB语言编写。这些代码远未完善，可修改余地很大，仅供教学之用。

function gradf=grad(hfun,x)
%GRAD  数值法求函数在给定点处的导数值（一元函数）或梯度（多元函数）
%   gradf = grad(hfun,x0)  hfun是函数句柄或内联函数，x0是一定点或一批点（按列）；返回
%   值gradf是函数在该点处的导数或梯度。
%   要求函数能对成批的点求函数值。比如：feval(hfun,X)返回一个与X同列的行向量，对应于以
%   X每一列作为函数自变量而求得的函数值。
%
%   Reference: 《最优化计算原理与算法程序设计》， 粟塔山等编著， 国防科技大学出版社
%  （湖南）， 2001.
%
%   $Author: WBC $    $Date: 2003/10/25 $

n = length(x);
h=1e-3;     % 数值法求梯度的步长
w1=zeros(n,1);%h/2
w2=w1; %-h/2
w3=w1; %h
w4=w1; %-h
for i=1:n
    x(i)=x(i)+h/2;
    w1(i)=hfun(x);
    x(i)=x(i)-h;
    w2(i)=hfun(x);
    x(i)=x(i)-h/2;
    w4(i)=hfun(x);
    x(i)=x(i)+2*h;
    w3(i)=hfun(x);
    x(i)=x(i)-h;
end
gradf=(8*(w1-w2)-(w3-w4))/(6*h);

function [fv,x,lambda,exitflag]=simplelinesearch(hf,rho,l,u,lambda0,fv0,x0,g,d)
%LINESEARCH 简单线搜索
%输入参数：
%        hf -- 函数句柄，目标函数
%       rho -- 实标量，简单线搜索的参数，介于0到0.5之间的数
%         l -- 实标量，简单线搜索的参数，介于0到1之间的数
%         u -- 实标量，简单线搜索的参数，介于0到1之间的数，满足u>l
%    alpha0 -- 实标量，步长的初始值
%       fv0 -- 实标量，一维搜索的初始目标函数值，即 hf(x0+alpha_0*d)
%        x0 -- 实列向量，当前点
%         g -- 实列向量，函数hf在当前点处的梯度
%         d -- 实列向量，函数hf在当前点处的搜索方向
%输出参数：
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
%          x -- 实列向量，下一个点
%     lambda -- 实标量，可接受步长
%   exitflag -- 整型标量，等于0表示线搜索成功，等于-1表示线搜索失败（内部迭代次数大于iterMax）

%参考文献：倪勤，最优化方法与程序设计，科学出版社
% Date: 2009/12/20

lambda = lambda0;
x = x0;
i = 0;
imax = 30; % 最大迭代次数，用户可以修改
% 主循环
gd = dot(g, d);
while i <= imax
    fv = hf(x + lambda*d);
    i = i + 1;
    if fv < fv0 + lambda * rho * gd;
        x = x + lambda*d;
        exitflag = 0;
        return;
    end
    lambda_bar = -gd*lambda^2*0.5/(fv-fv0-lambda*gd);
    lambda_bar = max(lambda_bar, l*lambda);
    lambda = min(lambda_bar, u*lambda);
end
exitflag=0;
if i >= imax && fv >= fv0
    fv = fv0;
    x = x0;
    lambda = 0;
    exitflag = -1;
end

function [fv, x, exitflag] = bfgs(hf, x0, epsi)
%BFGS 无约束问题的BFGS算法
%输入参数：
%        hf -- 函数句柄，目标函数
%        x0 -- 实列向量，初始点
%      epsi -- 实标量，终止误差
%输出参数：
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
%          x -- 实列向量，下一个点
%   exitflag -- 整型标量，等于0表示成功，等于-1表示失败（迭代次数大于iter_max）

%参考文献：倪勤，最优化方法与程序设计，科学出版社
% Date: 2009/12/20

%
%初始化
%
k = 0;
rho = 0.01;
l = 0.15;
u = 0.85;
x =x0;
fv = hf(x);
n = length(x);
H = eye(n);
iter_max = 100;
%
%检查终止条件
%
g = grad(hf, x);
while norm(g) > epsi && k<= iter_max
    d = -H*g;
    lambda0 = 1.0;
    %
    %做线搜索，如果成功，则返回更新当前点
    %
    [fv, x, lambda, exitflag] = simplelinesearch(hf, rho, l, u, lambda0, fv, x, g, d);
    if exitflag == -1 %重开始
        H = eye(n);
    else
        %
        %更新H
        %
        g_old = g;
        g = grad(hf, x);
        p = lambda*d;
        q = g- g_old;
        Hq = H*q;
        pq = p'*q;
        qHq = q'*Hq;
        v = sqrt(qHq) * (p/pq-Hq/qHq);
        H = H + p*p'/pq - Hq*Hq'/qHq + v*v';
    end
    k = k+1;
end
exitflag = 0;
if k > iter_max
    exitflag = -1;
end

function [fv, x, exitflag] = phr(hf, cf, x0)
%PHR  一般约束优化问题的PHR算法
%输入参数：
%         hf -- 函数句柄，目标函数
%         cf -- 函数句柄，约束条件，包含等式约束和不等式约束
%         x0 -- 实列向量，初始点
%输出参数：
%         fv --  实标量，一维搜索完成后的目标函数值，即 hf(x0+alpha*d)
%          x -- 实列向量，下一个点
%   exitflag -- 整型标量，等于0表示成功，等于-1表示失败（迭代次数大于iter_max）

%参考文献：倪勤，最优化方法与程序设计，科学出版社
% Date: 2009/12/20

%
%初始化
%
epsi = 1.0e-4;
k = 0;
sigma = 0.8;
c = 1.5;
theta = 0.8;
x = x0;
[ce, ci] = cf(x);
l = length(ce);
li = length(ci);
lambda = ones(l+li, 1) * 0.1;
iter_max = 100;
phi = 0;
if l
    phi = phi + ce'*ce;
end
if li
    phi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2);
end
while phi > epsi && k <= iter_max
    %
    %求解子问题
    %
    hmf = @(x) mfun(x, hf, cf,lambda, sigma);
    [fv, x] = bfgs(hmf, x, epsi);
    [ce, ci] = cf(x);
    phi_old = phi;
    phi = 0;
    if l
        phi = phi + ce'*ce;
    end
    if li
        phi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2);
    end
    if phi > epsi
        %
        %更新罚因子
        %
        if k >= 2 && phi/phi_old > theta
            sigma = c * sigma;
        end
        %
        %更新乘子
        %
        if l
            lambda(1:l) = lambda(1:l) - sigma*ce;
        end
        if li
            lambda(l+1:end) = max(0, lambda(l+1:end) - sigma*ci);
        end
    end
    k = k+1;
end
exitflag = 0;
fv = hf(x);
if k > iter_max
    exitflag = -1;
end

%
%乘子罚函数
%
function fv = mfun(x, hf, cf, lambda, sigma)
[ce, ci] = cf(x);
l = length(ce);
li = length(ci);
fv = 0;
fv = fv + hf(x);
if l
    fv = fv - lambda(1:l)'*ce + 0.5*sigma*ce'*ce;
end
if li
    fv = fv + 0.5/sigma*sum(max(0,lambda(l+1:end) - sigma*ci).^2 - lambda(l+1:end).^2);
end

 

 

这里是演示代码：

%% bfgs演示
%教材P328.1-3
hf1 = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; %banana函数
hf2 = @(x) (6 + x(1) + x(2)).^2 + (2 - 3*x(1) - 3*x(2) - x(1)*x(2)).^2;
hf3 = @(x) x(1).^2 - 2*x(1)*x(2) + 4*x(2).^2 + x(1) - 3*x(2);
[fv,x,exitflag]=bfgs(hf1,[0;0],0.001);
fv
x
exitflag
[fv,x,exitflag]=bfgs(hf2,[4;6],0.001);
fv
x
exitflag
[fv,x,exitflag]=bfgs(hf3,[1;1],0.001);
fv
x
exitflag
%% phr算法
%教材P414.5
hf1 = @(x) x(1).^2 + x(2).^2;
cf1 = @(x) deal([], x(1) - 1);
hf2 = @(x) x(1) + (x(2) + 1).^2/3;
cf2 = @(x) deal([],[x(1); x(2) - 1]);
%教材P392.2
hf3 = @(x) x(1).^2 + x(1).*x(2) + 2*x(2).^2 - 6*x(1) - 2*x(2) - 12*x(3);
cf3 = @(x) deal(x(1)+x(2)+x(3)-2,...
    [x(1) - 2*x(2) + 3; x(1); x(2); x(3)]);
%其它例子
hf4 = @(x) 6*x(2)*x(5) + 7*x(1)*x(3) + 3*x(2)^2;
cf4 = @(x) deal([3*x(2)^2*x(5) + 3*x(1)^2*x(3) - 20.875;
    x(1) - 0.3*x(2)],...
    [-x(1) + 0.2*x(2)*x(5) + 71
     -0.9*x(3) + x(4)^2 + 67
     x(3)
     x(5) - 1
     -x(3) + 20
     x(4) - 0.1*x(5)
     -x(4) + 0.5*x(5)
     x(3) - 0.9*x(5)
     ]);
hf5 = @(x) exp(x(1)) * (4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) + 2*x(2) + 1);
cf5 = @(x) deal([], [x(1) + x(2) - x(1)*x(2) - 1.5; x(1)*x(2) + 10]);

[fv, x, exitflag] = phr(hf1, cf1, [3;2]);
fv
x
exitflag
[fv, x, exitflag] = phr(hf2, cf2, [3;2]);
fv
x
exitflag
[fv, x, exitflag] = phr(hf3, cf3, [1;1;0]);
fv
x
exitflag
[fv, x, exitflag] = phr(hf4, cf4, [1; 4; 5; 2; 5]);
fv
x
exitflag
[fv, x, exitflag] = phr(hf5, cf5, [-1; 1]);
fv
x
exitflag