外罰函數主要用於對於等式約束問題的求解，內點法主要是對於不等式問題的求解，一般問題中包含等式約束以及不等式約束，故需要使用乘子法解決問題。
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=-Bk\gk;%解方程組，計算搜索方向
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)
輸出結果如下：
