乘子法matlab

运用乘子法求解问题

比如这个

通过构造几个函数

然后带入已知量

clc,clear;
clear;
clc;
x0=[0,0,0]';
[x,mu,lam,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0)
函数：
function [x,mu,lam,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)
maxk=1000;
sigma=2.0;
theta=0.8; 
eta=2.0; 
k=0; 
ink=0;
epsilon=1e-5;
x=x0; 
he=feval(hf,x); 
gi=feval(gf,x);
n=length(x); 
l=length(he); 
m=length(gi);
mu=0.1*ones(l,1); 
lam=0.1*ones(m,1);
betak=10; 
betaold=10;
while(betak>epsilon && k<maxk)
    [ik,x,val]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lam,sigma);
    ink=ink+ik;
    he=feval(hf,x); 
    gi=feval(gf,x);
    betak=sqrt(norm(he,2)^2+norm(min(gi,lam/sigma),2)^2);
    if betak>epsilon
        mu=mu-sigma*he;
        lam=max(0.0,lam-sigma*gi);
        if(k>=2 && betak>theta*betaold)
            sigma=eta*sigma;
        end
    end
    k=k+1;
    betaold=betak;
    x0=x;
end
f=feval(fun,x);
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.beta=betak;
 

函数：
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lam,sigma)
f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);
l=length(he); m=length(gi);
psi=f; s1=0.0;
for i=1:l
    psi=psi-he(i)*mu(i);
    s1=s1+he(i)^2;
end
psi=psi+0.5*sigma*s1;
s2=0.0;
for i=1:m
    s3=max(0.0,lam(i)-sigma*gi(i));
    s2=s2+s3^2-lam(i)^2;
end
psi=psi+s2/(2.0*sigma);
 

函数：
function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lam,sigma)
dpsi=feval(dfun,x);
he=feval(hf,x); 
gi=feval(gf,x);
dhe=feval(dhf,x); 
dgi=feval(dgf,x);
l=length(he); 
m=length(gi);
for(i=1:l)
    dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);
end
for(i=1:m)
    dpsi=dpsi+(sigma*gi(i)-lam(i))*dgi(:,i);
end

函数：
function[k,x,val]=bfgs(fun,gfun,x0,varargin)
N=1000;
epsilon=1.e-5;
beta=0.55;
sigma=0.4;
n=length(x0);
Bk=eye(n);
k=0;
while(k<N)
    gk=feval(gfun,x0,varargin{:});
    if(norm(gk)<epsilon)
        break;
    end
    dk=-Bk\gk;
    m=0;
    mk=0;
    while(m<20)
        newf=feval(fun,x0+beta^m*dk,varargin{:});
        oldf=feval(fun,x0,varargin{:});
        if(newf<=oldf+sigma*beta^m*gk'*dk)
            mk=m;
            break;
        end
        m=m+1;
    end
    x=x0+beta^mk*dk;
    sk=x-x0;
    yk=feval(gfun,x,varargin{:})-gk;
    if(yk'*sk>0)
        Bk=Bk-(Bk*(sk*sk')*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);
    end
    k=k+1;
    x0=x;
end
val=feval(fun,x0,varargin{:});
 

带入一些已知条件
function f=f1(x)
%要求的最小值函数

    f=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(3);
end
function he=h1(x)
%等式约束

    he=[8*x(1)+14*x(2)+7*x(3)-56,x(1)^2+x(2)^2+x(3)^2-25]';
 
end
function gi=g1(x)
%不等式约束

gi=zeros(3,1);
gi(1)=x(1);
gi(2)=x(2);
gi(3)=x(3);
 
end
function g=df1(x)
%要求的函数求导哦

    g=[-2*x(1)-x(2)-x(3),-4*x(2)-x(1),-2*x(3)-x(1)]';
end
function dhe=dh1(x)
%等式约束求导

    dhe=[8,14,7;2*x(1),2*x(2),2*x(3)]';
end
function dgi=dg1(x)
%不等式约束求导

dgi=[1 0 0;0 1 0;0 0 1];
end

 

就是这样