罚函数罚与乘子法 (1) 2

罚函数法

罚函数法是能够处理一般的约束优化问题：min ()()0,1,2,()0,1,2,,i i

f x h x i k

g x j m ??

==??≥=?的一类方

法。其基本思想是将约束优化问题卑微无约束问题来求解。罚函数是由目标函数和约束函数的某种组合得到的函数，对于等式约束的优化问题

min ()()0,1,2,

i f x h x i k

??

==?，可以定义如下的罚函数：

21

()()()k

i i F x f x c h x ==+∑

将约束优化问题转化为无约束优化问题；对于不等式约束的优化问题

min ()()0,1,2,

,i f x g x j m

??

≥=?

可以定义如下的罚函数：

1

1

()()()

m

j j F x f x C g x ==+∑

对于同时存在等式约束和不等式约束的优化问题，可以去上面两个罚函数的组合。当然罚函数还有其他的取法，但是构造罚函数的思想都是一样的，即使得在可行点罚函数等于原来的目标函数值，在不可行点罚函数等于一个很大的数。

外点罚函数法 1.算法原理

外点罚函数法是通过一系列罚因子{}i c ，求罚函数的极小值来逼近原约束问题的最有点。之所以称为外点罚函数法，是因为它是从可行域外部向约束边界逐步靠拢的。 2，。算法步骤

用外点罚函数法求解线性约束问题min ()

f x Ax b ??=?的算法过程如下：

1，给定初始点(0)x ，罚参数列{}i c 及精度0ε>，置1k =； 2，构造罚函数2

()()F x f x c Ax b =+-；

3，用某种无约束非线性规划，以(1)k x -为初始点求解min ()F x ；

4，设最优解为()k x ，若()k x 满足某种终止条件，则停止迭代输出()k x ，否则

令1k k =+，转2；

罚参数列{}i c 的选法：通常先选定一个初始常数1c 和一个比例系数2ρ≥，则其余的可表示为11i i c c ρ-=。终止条件可采用()S x ε≤，其中2

()S x c Ax b =-。 3 算法的MATLAB 实现

function [x,minf] = minPF(f,x0,A,b,c1,p,var,eps) % 目标函数：f ； % 初始点： x0; % 约束矩阵： A ； % 约束右端向量：b ； % 罚参数的初始常数： c1; % 罚参数的比例系数： p ； % 自变量向量 Var ； % 精度： eps ；

% 目标函数取最小值时自变量值： x; % 目标函数的最小值：minf ； format long ; if nargin == 7 eps = 1.0e-4; end k = 0; FE = 0;

for i=1:length(b)

FE = FE + (var*transpose(A(1,:)) - b(i))^2; end

x1 = transpose(x0); x2 = inf; while 1 M = c1*p; FF = M*FE; SumF = f + FF;

[x2,minf] = minNT(SumF,transpose(x1),var); if norm(x2 - x1)<=eps x = x2; break ; else c1 = M; x1 = x2; end end

minf = Funval(f,var,x); format short ;

4.算法举例

用外点罚函数法求下面的优化问题：22min (,)0.50.25 ,.. 1f t s t s s t t s =++=其中取10.05,2c ρ==，初始点取0(0,0)x =。

解：由题已知[1 1];1A b ==，则每一步迭代需求解的罚函数为：

222(,)0.50.250.05*2(1)i f t s t s t s =+++-

在MATLAB 命令窗口中输入： >> f = 0.5*t^2+s^2/4; >> A=[1 1];b=1; >> c1=0.05;p=2;

>> [x,minf]=minPF(f,[0 0],A,b,c1,p,[t s]) 所得结果为： >> x=0.3333 0.6666 minf = 0.1666

对于一般的等式约束问题min ()

()0

f x h x ??=?也可以用外点罚函数法解决：

function [x,minf] = minGeneralPF(f,x0,h,c1,p,var,eps) format long ; if nargin == 6 eps = 1.0e-4; end k = 0; FE = 0;

for i=1:length(h) FE = FE + (h(i))^2; end

x1 = transpose(x0); x2 = inf; while 1 M = c1*p; FF = M*FE; SumF = f + FF;

[x2,minf] = minNT(SumF,transpose(x1),var); if norm(x2 - x1)<=eps x = x2;

else c1 = M; x1 = x2; end end

minf = Funval(f,var,x); format short ;

例 用通用罚函数法求下面的优化问题：

min (,)2, .. 1, ,0f t s t s s t ts t s =+=> 其中取10.05,2c p ==，初始点取为

0(0,0)x =。

外点罚函数法除了可以求含等式约束的优化问题外，还可以求混合约束的优化问题，不过其罚函数形式要修正。

min ()()0,1,2,()0,1,2,,i i

f x h x i k

g x j m ??

==??≥=? 通过引入罚函数

[]2

21

1

() ()()min(0,())k m

i i i j F x f x M h x M g x ===++∑∑

可以转化为无约束优化问题。

function [x,minf] = minConPF(f,x0,g,h,c1,p,var,eps)

format long ; if nargin == 7 eps = 1.0e-4; end k = 0; FE = 0;

for i=1:length(h) FE = FE + (h(i))^2; end

x1 = transpose(x0); x2 = inf; while 1 M = c1*p;

gx = Funval(g,var,x1); gF = 0;

for i=1:length(h) if gx(i) < 0

gF = gF + M*(g(i)^2); %构造罚函数 end end

SumF = f + FF + gF;

[x2,minf] = minNT(SumF,transpose(x1),var); %用牛顿法求解无约束规划问题

if norm(x2 - x1)<=eps x = x2; break ; else c1 = M; x1 = x2; end end

minf = Funval(f,var,x); format short ;

内点罚函数法 1.算法原理

内点罚函数法的所有迭代过程均在可行域内进行，它每次迭代得到的点都是可行点，内点罚函数法的基本思想是用可行域内的点列来逼近最优解。 2，。算法步骤

用内点罚函数法求解线性约束问题min ()

()0,1,2,,

i f x g x j m ??

≥=?的算法过程如

下：

1，给定初始点(0)x ，罚因子0μ>，缩小系数1v >及精度0ε>，置1k =； 2，构造增广目标函数()()()F x f x B x μ=+，其中11

()()

m

i i B x g x ==∑

； 3，用某种无约束非线性规划，以(1)k x -为初始点求解min ()F x ；设最优解为()k x ，若()()k B x με<，则终止迭代输出()k x ，否则令v μμ=，1k k =+转2； 3 算法的MATLAB 实现

function [x,minf] = minNF(f,x0,g,u,v,var,eps)

format long; if nargin == 6 eps = 1.0e-4; end k = 0; FE = 0;

for i=1:length(g)

FE = FE + 1/g(i); %构造罚函数 end

x1 = transpose(x0); x2 = inf;

while 1

FF = u*FE; SumF = f + FF;

[x2,minf] = minNT(SumF,transpose(x1),var); %用牛顿法求解无约束规划 Bx = Funval(FE,var,x2);

if u*Bx < eps %精度判断 if norm(x2 - x1)<=eps x = x2; break; else

u = v*u; %参数修正 x1 = x2; end else

if norm(x2 - x1)<=eps x = x2; break; else

u = v*u; %参数修正 x1 = x2; end end end

minf = Funval(f,var,x); format short;

4 算法举例

2211

min (,), .220

t s f t s t s s t t s +-≥?=+?

--≥? 其中取8,0.5v μ==，初始点取为

(,)(3,1)t s =。

>> syms t s; >> f = t^2+s^2;

>> g=[t+s-1;2*t-s-2];

>> [x, minf]=minNF(f, [3 1], g, 8, 0.5,[t s]) x=1.0002 0.0001 minf= 1.0004

乘子法

1.算法原理

外点罚函数法和内点罚函数法的均要求罚因子趋于无穷才能得到目标函数的最优解，但是罚因子太大的话会引起计算的困难，为克服这一缺陷就产生了乘子法。

2，算法步骤

用乘子法求解约束优化问题min ()

()0,1,2,,i f x h x j m

??

==?的算法过程如下：

1，给定初始点(0)x ，初始乘子向量12(,,,)m v v v v ，罚因子0M >及放大系

数1α>，及精度0ε>，参数(0,1)γ∈，置1k =；

2，构造目标函数21

1

()()()()2

m m

i

i i

i i M

F x f x h x v h x ===+

-∑∑；

3，用某种无约束非线性规划，以(1)k x -为初始点求解min ()F x ；设最优解为()k x ，

4，若()()k h x ε<，则终止迭代输出()k x ，否则转2； 5，若

()(1)

()()

k k h x h x

γ-≥，则令M M α=，否则保持M 不变，转6

6，令()(),1,2,,,k i i i v v Mh x i m =-= 置1k k =+，转2，

3 算法的MATLAB 实现

function [x,minf] = minFactor(f,x0,h,v,M,alpha,gama,var,eps) format long ; if nargin == 8

eps = 1.0e-4; end FE = 0;

for i=1:length(h) FE = h(i)^2; end

x1 = transpose(x0); x2 = inf; while 1

FF = M*FE/2; Fh = v*h;

SumF = f + FF - Fh;

[x2,minf] = minNT(SumF,transpose(x1),var);

Hx2 = Funval(h,var,x2); Hx1 = Funval(h,var,x1); if norm(Hx2) < eps x = x2; break ; else

if Hx2/Hx1 >= gama M = alpha*M; x1 = x2; else

v = v - M*transpose(Hx2); x1 = x2; end end end

minf = Funval(f,var,x); format short ;

4 算法举例

22min (,), . 21f t s t s s t t s =++= 其中取8,3,0.7,(1),M v αγ====，初始点取为0(0,0)x =。

>> syms t s; >> f = t^2+s^2;

>> h=[t+s-1;2*t-s-2];

>> [x, minf]=minFactor(f, [0 0], h, 2,3, 0.7,[t s])

x=0.2000 0.4000 minf= 0.2000