约束优化算法拉格朗日乘子法

拉格朗日乘子法约束优化问题的标准形式为:min (),..()0,1,2,...,()0,1,2,...,ni j f x x R s t g x i m h x j l∈≤===,,:n i j f g h R R →其中约束优化算法的基本思想就是:通过引入效用函数的方法将约束优化问题转换为无约束问题,再利用优化迭代过程不断地更新效用函数,以使得算法收敛。

1. 罚函数法罚函数法(内点法)的主思想就是:在可行域的边界上筑起一道很高的“围墙”,当迭代点靠近边界时,目标函数陡然增大,以示惩罚,阻止迭代点穿越边界,这样就可以将最优解“挡”在可行域之内了。

它只适用于不等式约束:min (),..0,1,2,...,ni f x x R s tg i m ∈≤=它的可行域为: {|()0,1,2,...,}n i D x R g x i m =∈≤=对上述约束问题,其其可行域的内点可行集0D ≠∅的情况下,引入效用函数:min (,)()()B x r f x rB x =+%、 其中11()()m i i B x g x ==-∑%或1()|ln(())|m i i B x g x ==-∑% 算法的具体步骤如下:给定控制误差0ε>,惩罚因子的缩小系数01c <<。

步骤1:令1k =,选定初始点(0)0x D ∈,给定10r >(一般取10)。

步骤2:以()k x 为初始点,求解无约束min (,)()()k B x r f x r B x =+% 其中11()()m i i B x g x ==-∑%或1()|ln(())|m i i B x g x ==-∑%,得最优解()()k k x x r = 步骤3:若()()k k r B x ε<%,则()k x 为其近似最优解,停;否则,令,1k k r cr k k ==+,转步骤2、2. 拉格朗日乘子法(1)PH 算法:(约数为等式的情况引入)效用函数为()()min (,,)()()()()k k T T k k M x u f x u h x h x h x σσ=++判断函数为()()k k h x φ=当()()k k x φφε=<时迭代停止。

步骤1:选定初始点(0)x ,初始拉格朗日乘子向量(1)u ,初始罚因子1σ及其放大系数1c >,控制误差0ε>与常数(0,1)θ∈,令1k =。

步骤2:以(1)k x +为初始点,求解无约束问题:()()min (,,)()()()()k k T T k k M x u f x u h x h x h x σσ=++得到无约束问题最优解()k x步骤3:当()()k h xε<时,()k x 为所求的最优解,停;否则转步骤4、 步骤4:当()()()()/k k h xh x θ<时,转步骤5;否则令k 1k c σσ+=,转步骤5、 步骤5:令(1)()()(),1k k k k u u h x k k σ+=+=+,转步骤1。

(2) PHR 算法(一般约束形式的松弛变量法与指数形式法)松弛变量法: (){}12222111(,,)()max 0,()2()()2i i i l lj j jj j M u v f x u g x u v h x h x ρρρρρ===++-⎡⎤⎣⎦++∑∑∑乘子的修正公式为:(1)()()(1)()()(),1,...,max 0,(),1,...,k k k j j j k k k i i i v v h x j luu g x i m ρρ++=+=⎡⎤=+=⎣⎦判断函数为: 1/22()2()()11()max (),k l m k k i k j i j i u h x g x φρ==⎧⎫⎛⎫⎪⎪=+-⎨⎬ ⎪⎝⎭⎪⎪⎩⎭∑∑ 当()()k k x φφε=<时迭代停止。

3. 乘子法MATLAB 程序及其作用Al main函数3、1 _3、1、1程序(1):乘子法效用函数程序函数功能:将约束优化问题,根据效用函数方法,将其转变成无约束问题。

function f=AL_obj(x)%拉格朗日增广函数%N_equ 等式约束个数%N_inequ 不等式约束个数global r_al pena N_equ N_inequ;%全局变量h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式约束部分for i=1:N_equh_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i)、^2;end%不等式约束部分for i=1:N_inequh_inequ=h_inequ+(0、5/pena)*(max(0,(r_al(i)+pena*g(i)))、^2-r_al(i)、^2); end%拉格朗日增广函数值f=obj(x)+h_equ+h_inequ;3.1.2程序(2):判断函数函数功能:判断就是否符合约束条件%% the compare function is the stop conditionfunction f=compare(x)global r_al pena N_equ N_inequ;h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式部分for i=1:N_equh_equ=h_equ+h(i)、^2;end%不等式部分for i=1:N_inequh_inequ=h_inequ+(max(-g(i),r_al(i+N_equ)/pena))、^2;endf=sqrt(h_equ+h_inequ);3.1.3程序(3)AL算法主程序函数功能:对无约束的效用函数利用拟牛顿算法求解其最优解,更新乘子。

function [X,FVAL]=AL_main(x_al,r_al,N_equ,N_inequ)%本程序为拉格朗日乘子算法示例算法%函数输入:% x_al:初始迭代点% r_al:初始拉格朗日乘子% N-equ:等式约束个数% N_inequ:不等式约束个数%函数输出% X:最优函数点% FVAL:最优函数值%============================程序开始================================ global r_al pena N_equ N_inequ; %参数(全局变量)pena=10; %惩罚系数c_scale=2; %乘法系数乘数cta=0、5; %下降标准系数e_al=0、005; %误差控制范围max_itera=25;out_itera=1; %迭代次数%===========================算法迭代开始============================= while out_itera<max_iterax_al0=x_al;r_al0=r_al;%判断函数compareFlag=compare(x_al0);%无约束的拟牛顿法BFGS[X,FVAL]=fminunc(@AL_obj,x_al0);x_al=X; %得到新迭代点%判断停止条件if compare(x_al)<e_aldisp('we get the opt point');breakend%c判断函数下降度if compare(x_al)<cta*compareFlagpena=pena; %可以根据需要修改惩罚系数变量elsepena=min(1000,c_scale*pena); %%乘法系数最大1000disp('pena=2*pena');end%% 更新拉格朗日乘子[h,g]=constrains(x_al);for i=1:N_equ%%等式约束部分r_al(i)=r_al(i)+pena*h(i);endfor i=1:N_inequ%%不等式约束部分r_al(i+N_equ)=max(0,(r_al(i+N_equ)+pena*g(i)));endout_itera=out_itera+1;end%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++ disp('!!!!!!!!!!!!!!!!!!!the iteration over!!!!!!!!!!!!!!!!!!!!!!!!!!');disp('the value of the obj function');obj(x_al)disp('the value of constrains');compare(x_al)disp('the opt point');X=x_al;FVAL=obj(X);3.1.4 乘子法_AL main 函数使用方法（1） 定义目标函数及约束条件122331123222123min ()..1030f x x x x x x x s tx x x x x x =---++-=++-≤目标函数m 文件 function f=obj(x)f=-x(1)*x(2)-x(2)*x(3)-x(3)*x(1);约束函数m 文件222[,]()(1)(2)(3)1;(1)(2)(3)3;function h g constrains x h x x x g x x x ==++-=++-（2） _AL main 函数调用x_al=[1,1,1]; %初始迭代点r_al=[1,1]; %初始拉格朗日乘子N_equ=1; %等式约束个数 一个N_inequ=1; %不等式约束个数 一个[X,FVAL]=AL_main(x_al,r_al,N_equ,N_inequ)计算结果:we get the opt point!!!!!!!!!!!!!!!!!!!the iteration over!!!!!!!!!!!!!!!!!!!!!!!!!!the value of the obj functionans =-3、9871e+031the value of constrainsans =the opt pointX =1、0e+015 *3、7723 3、3985 3、7723 FVAL =-3、9871e+031。