外罚函数主要用于对于等式约束问题的求解，内点法主要是对于不等式问题的求解，一般问题中包含等式约束以及不等式约束，故需要使用乘子法解决问题。
1、 乘子法概述
（1）等式约束乘子法描述：
min f(x)
s.t. gi(x) =0
广义乘子法是拉格朗日乘子法与罚函数法的结合，构造增广函数：
φ (x,λ,σ)=f(x)+λTg(x)+1/2σgT(x)g(x)
在罚函数的基础上增加了乘子项，首先在σ足够大的基础上，获得ϕ的极小值，然后在调整λ获得原问题的最优解。
（2）包含等式约束以及不等式约束问题描述：
min f(x)
s.t. hi(x) =0,i=1,...,l
gi(x)≥0,i=1,...m
其基本思想是：先引进辅助变量把不等式约束化为等式约束，然后利用最优性条件消去辅助变量，主要是通过构造增广拉格朗日函数，进行外迭代与内迭代综合，带入乘子迭代公式，进而得出得出，故针对上述一般问题构造拉格朗日函数为：
4、其代码实现为
function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)
%功能：用乘子法解一般约束问题：min f(x),s.t. h(x)=0.g(x)>=0
%输入：x0是初始点，fun，dfun分别是目标函数及其梯度；
%hf，dhf分别是等式约束（向量）函数及其jacobi矩阵的转置；
%gf，dgf分别是不等式约束（向量）函数及其jacobi矩阵的转置；
%输出：x是近似最优点，mu，lambda分别是相应于等式约束和不等式
% 等式约束的乘子向量；output是结构变量，输出近似极小值f，迭代次数，内迭代次数等
%%%%%%c初始化相关参数%%%%%%%%%%%
maxk=500; %最大迭代次数
sigma=2.0; %罚因子
eta=2.0; theta=0.8; %PHR算法中的实参数
k=0; ink=0; %k,ink分别是外迭代和内迭代次数
epsilon=1e-5;%终止误差值
x=x0;he=feval(hf,x);gi=feval(gf,x);%he=feval(hf,x)=hf(x)
n=length(x);l=length(he);m=length(gi);
%选取乘子向量的初始值
mu=0.1*ones(1,1);lambda=0.1*ones(m,1);%ones为生成m*n的全1矩阵
btak=10; btaold=10; %用来检验终止条件的两个值
while (btak>epsilon & k<maxk)
%%%%%%c先求解无约束问题%%%%%%%%%%%
%调用BFGS算法程序求解无约束子问题
[x,v,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma);%%其中x为最优点，val为最优值，ik为迭代次数
ink=ink+ik;
he=feval(hf,x);gi=feval(gf,x);
%%%%%%%%%%计算btak%%%%%%%%%%%
btak=0.0;
for(i=1:l),btak=btak+he(i)^2; end
for(i=1:m)
temp=min(gi(i),lambda(i)/sigma);
btak=btak+temp^2;
end
btak=sqrt(btak);
if btak>epsilon
%%%%%%%%%%%更新罚参数%%%%%%%%%%%
if(k>=2 & btak>theta*btaold)
sigma=eta*sigma;
end
%%%%%%%%%%%更新乘子向量%%%%%%%%%%%%
for(i=1:l),mu(i)=mu(i)-sigma*he(i);end
for(i=1:m)
%lambda(i)=max(0.0,lambda(i)-sigma*gi(i));
lambda(i)=max(0.0,lambda(i)-gi(i));
end
end
%%%%%%%%%%%迭代%%%%%%%%%%%%
k=k+1;
btaold=btak;
x0=x;
end
f=feval(fun,x);
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.bta=btak;
BFGS算法部分：
function [x,val,k]=bfgs(fun,gfun,x0,varargin)
%功能：用BFGS算法求解无约束问题：minf（x）
%输入：x0是初始点，fun，gfun分别是目标函数及其梯度
%varargin是输入的可变参数变量，简单调用bfgs时可以忽略，其他程序调用则尤为重要
%输出：x为最优点，val为最优值，k时迭代次数
maxk=500;%给出最大迭代次数
rho=0.55;sigma=0.4;epsilon=1e-5;
k=0;n=length(x0);
Bk=eye(n);%Bk=feval('Hess',x0)
while(k<maxk)
gk=feval(gfun,x0,varargin{:});%计算梯度
if(norm(gk)<epsilon),break;end%检验终止准则
dk=-Bkgk;%解方程组，计算搜索方向
m=0;mk=0;
while(m<20)%搜索求步长
newf=feval(fun,x0+rho^m*dk,varargin{:});
oldf=feval(fun,x0,varargin{:});
if(newf<oldf+sigma*rho^m*gk'*dk)
mk=m;break;
end
m=m+1;
end
%bfgs校正
x=x0+rho^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=(x(1)-2.0)^2+(x(2)-1.0)^2;
%等式约束条件
function he=h1(x)
he=x(1)-2.0*x(2)+1.0;
%不等式约束条件
function gi=g1(x)
gi=-0.25*x(1)^2-x(2)^2+1;
%目标函数的梯度文件
function g=df1(x)
g=[2.0*(x(1)-2.0),2.0*(x(2)-1.0)]';
%等式函数的Jacobi（转置）矩阵文件
function dhe=dh1(x)
dhe=[1.0,-2.0]';
%不等式约束函数的Jacobi矩阵（转置矩阵）
function dgi=dg1(x)
dgi=[-0.5*x(1),-2.0*x(2)]';
命令行指令为：
x0=[3,3]'
[x,mu,lambda,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0)
输出结果如下：
