最优化理论与方法

内点法基本原理

摘要：内点法是求解含不等式约束最优化问题的一种十分有效的算法。内点法通过构造障碍函数，求解一系列只含等式约束最优化问题，逐步得到原问题的最优解，具有找初始点容易、线性收敛、迭代次数少等特点。本文主要介绍了内点法的基本原理，障碍方法的一般步骤并分析了该方法的优缺点，进行了算例实践。

关键词：内点法；障碍方法；Newton法

The Theory of Interior Point Method

Abstract: Interior point method is a very effective algorithm for solving optimization problems with inequality constrained. Interior point method is constructed to solve a series of optimization problems with equality constraints, and the optimal solution of the original problem is obtained, which has the characteristics of finding the initial point easier, linear convergence, less iteration number and so on.

This paper mainly introduces the theory of interior point method, the general steps of barrier method and analyzing the advantages and disadvantages of the method.

Key words: interior point method; barrier method;Newton method

1 引言

内点法是由John von Neumann 利用戈尔丹的线性齐次系统提出的，后被Narendra Karmarkar 于1984年推广应用到线性规划，即Karmarkar 算法。内点法是凸优化算法递阶结构中新层次的方法，其思想是对于线性等式和不等式约束的优化问题，通过将其简化成一系列线性等式约束问题求解[1]。主要是通过构造对数障碍函数将不等式约束与原目标函数结合，将原问题转化为新的等式约束下的优化问题，之后再运用Newton 方法进行求解[2]。而罚函数法是将不等式约束与等式约束全部进行加权处理后与原目标函数结合，将原问题转化为新的无约束优化问题，通过求解该新的无约束优化问题，间接得到原约束优化问题的最优解[6]。两者在算法思想上有本质的不同，但是当原问题中只有不等式约束时，两者求解是相同的。下面主要介绍内点法的基本原理及算法思想，并对障碍法与原对偶内点法的优缺点进行了探究与分析，并给出了一个实际算例来加以说明。

2 内点法原理

2.1 对数障碍函数和中心路径

含有不等式约束的凸优化问题的标准形式如下：

b

Ax m i x f x f i ==≤,,1,0)(subject to )

(minimize 0 （1）

其中R R f f n m →:,,0 是二次可微的凸函数，n

p R A ?∈，n p rankA <=。并且设最优的点*

x ，最优值

为)(0**=x f p 。满足如下KKT 条件：

m

i x f A x f x f m

i m i x f b Ax i i

m

i T i i i i 2,1,

0)(0

)()(2,10

2,1,

0)(,10==?=+?+?=>=≤=**=*******∑λυλλ（2）

通过定义对数障碍函数及加入正参数t ：

,},1,0)(|{)()),(log()(1

∑==<∈=--=m

i i n i m i x f R x x dom x f x φφ

将原问题（1）近似转化为

b

Ax x x f t =+subject to

)()(minimize 10φ（3）

随着参数t 的不断增加问题（3）的近似精度也会不断改进。而后我们可以通过让参数t 逐渐增加，运用Newton 方法来求解一系列形如（3）的优化问题，从而得到原问题（1）的最优解。 为了简化符号，考虑（3）的等价问题：

b

Ax x x tf =+subject to )

()(minimize 0φ（4）

两则的最优解集相同，最优值差一个常参数t 。从开始假定问题（4）能够用Newton 方法求解，并特别假定对任何t>0都存在唯一解)(t x *。我们称)(t x *为中心点，而这些点的集合定义为问题（1）的中心路径。并且所有中心路劲上的点都由以下充要条件所界定：)(t x *是严格可行，即满足：

,2,1,0))((,)(m i t x f b t Ax i =<=**

并且存在p

?R ∈υ

使 υ

φ?))(())((00T A t x t x f t +?+?=** 成立。再定义

,/?)( ,2,1,))

((1)(t t m i t x tf t i i υυλ==-

=*

*

* 我们可以将中心路径的条件解释为KKT 最优性条件（2）的连续变形。即x 等于)(t x *的充要条件是存在

υλ，满足：

m

i t x f A x f x f m

i m i x f b Ax i i m

i T

i i i i 2,1,

1/)(0

)()(2,102,1,

0)(,1

0==?=+?+?=>=≤=∑=λυλλ（5）

中心路径条件（5）与KKT 条件（2）的唯一不同在于互补松弛条件0)(-=x f i i λ被条件t x f i i 1/)(-=λ所替换。特别是，对于很大的t ，)(t x *

和对应的对偶解)()(t t *

*

υλ，“几乎”满足问题（1）的KKT 最优性条件。并且在实际计算过程中求解0)(-=x f i i λ是相当有难度的，而通过将其转化为t x f i i 1/)(-=λ就相对容易求解。

2.2 障碍方法

障碍方法是对无约束极小化方法进行简单的扩充，通过顺序求解一系列无约束（或线性约束）对的极小化问题，每次用所获得的最新的点作为求解下一个问题的初始点，也就是说，通过逐渐增加t 的值得到相应的)(t x *

，直到满足所需要的精度。 障碍算法的步骤如下：

给定严格可行点，

1,0:,0

>>=μt t x 误差阈值0>ε。

重复进行：

1.中心点步骤。

从x 开始，在b Ax =的约束下极小化φ+0tf ，最终确定)(t x *。 2.改进。)(:t x x *=。

3.停止准则。如果对偶间隙小于ε。

4.增加t 。t t μ=:。

每步迭代中我们都从上次获得的中心点开始计算当前中心点)(t x *，然后通过乘比因子1>μ增加t 。在中心点步骤中我们主要采用Newton 方法来求解。我们将中心点步骤中的Newton 迭代或步骤称为内部迭代。每次内部迭代我们都可以得到原问题可行解；但是仅在外部迭代结束时我们才有对偶可行解。 2.3 障碍方法优缺点分析 （1）对初始点选择要求不高

对于初始点的选择可以通过求解下面优化问题得到，

b

Ax m i s x f s

i ==≤ 2,1,)(subject to minimize （6）

并且根据上述优化问题的最优值*

p 的符号可以区分三种情况。当0<*p 则原不等式约束有严格可行解；当0>*p 则原不等式约束不可行；当0=*

p 则原不等式约束可行，但是没有严格可行解。

并且只需要选择的初始点严格满足不等式约束，在接下去的迭代中只要中心步骤的某个迭代点选取了完整的Newton 步径，那么之后的所有迭代点都是原问题可行点。 （2）中心点的精度不需要十分精确

在求解过程中不需要对)(t x *

进行精确计算，因为中心路径的作用仅是随着∞→t 将中心点引向原点的最优解，无其他意义；不精确的中心点计算方法同样能够产生收敛于最优点的点列)

(k x

。并且另一个方

面，计算φ+0tf 的十分精确的极小解只会增加少量的Newton 迭代，因此也可以假定在计算中心点步骤中产生的都是精确的中心点。 （3）μ值的选择不敏感

参数μ的选择需要同时兼顾所需要的内部迭代和外部迭代的次数。当μ值较小时可以期望每次外部迭代所需要进行较少次数的Newton 迭代，但是由于每次外部迭代只减少了较小的间隙，所需要的外部迭代次数会比较多。当μ较大时，相反的情况也会发生。但是经过实践表明，μ的选择对总的计算量并非特别敏感，取值从10到20或附近都能取得较好的效果，这位算法选择μ提供了极大的便利。 （4）随着问题维数的增加所需要的Newton 迭代次数的增量非常小

在具体算例中，可以发现应用障碍方法当问题的维数增加时，所需要的Newton 迭代次数非常缓慢地增加，几乎总是围绕数十次的数量变动。不过，进行每次Newton 迭代的计算量会随着问题规模的增加而同步增加。

（5）对偶间隙呈现近似线性收敛

对于大于10或者附近数值的μ值，障碍方法能够很好地在大约经过30次左右的Newton 迭代就可以把对偶间隙从2

10减少到-6

10。每次Newton 迭代问题的对偶间隙近似缩小为原来的μ1/。但是通过实践，我们也发现当μ大于10或附近的数值时，μ的选择几乎不影响所需要的总的Newton 迭代次数。为了改进障碍方法的线性收敛速度，又提出了原对偶内点法，该方法具有超线性收敛性质，并且能够有效处理可行但不严格可行的问题[5]。

3 算法实践

为了编程的方便，令t

r 1=，则0lim ,0k

1

=>>>>∞

→k

k r r r r 。给定精度等参数后按照上述算法可以得到如下程序框图：

从 点出发求解 ：

(0)1,,,r C

εε开始

在可行域内选取

(0)

X k=1

(1)k X -()*()min (,),X ()

k k X r r ?得

()()11

()()k k X r X r ε-**-≤()()1()()k k f X f X ε--≤ (1)()

k k r Cr += 1

k k =+

**()()

k X X r =**()()(())

k f X f X r =停止

图1：内点法程序框图

以下面简单的优化问题为例对内点进行进一步的说明。

1)(subject to )(minimize 112

2

21≤-=+=x x f x x x f （7）

首先通过构造内点障碍函数

()()

[][]22()1

211

(,)()ln ()ln (1)m

k k k u

u X r f X r

g X x

x r x ?==--=+---∑

用极值条件进行求解可得

()

111201

k r x x x ??=-=?-，2220x x ??==? 联立上式求得

()*()

1

112()2

k k r x r ±+=，*()2()0k x r =

由于约束条件的限制，可得无约束极值点为

()*()

112(),02T

k k r X r ??++=??????

，

当取0→k

r 时，可得最优解为

*[1,0]T X =，*()1f X =。

利用matlab 也得到了相同的结果，如附录程序所示[7]

。

而从如下图像中也可以直观看出当0→k

r 时求出的最优解约近似于原问题的最优解。

1

1.2

1.4

1.6

1.8

2

-2

-1

1

2

102030401

1.2

1.4

1.6

1.8

2

-2

-1

1

2

010

20

30

40

1

1.2

1.4

1.6

1.8

2

-2

-1

1

2

024681

1.2

1.4

1.6

1.8

2

-2

-1

1

2

02

4

6

8

图2：当r 分别取4、1、0.1、0.01时障碍函数图像

参考文献：

[1] Stephen Boyd, Lieven Vandenberghe.Convex Optimization[M].New York: Cambridge University Press, 2004,561-615.

[2] 刘红英, 夏勇, 周水生.数学规划基础[M].北京:北京航空航天大学出版社,2012,205-239.

[3] 刘明波, 王晓村. 内点法在求解电力系统优化问题中的应用综述[J]. 电网技术, 1999, 23(8):61-64.

[4] 汪超群, 韦化, 吴思缘,等. 七种最优潮流分解协调算法的性能比较[J]. 电力系统自动化, 2016, 40(6).

[5] 刘志鹏. 基于原对偶内点法的最优潮流研究[D]. 山东科技大学, 2009.

[6] 张菊亮, 章祥荪. 不等式约束最优化的非光滑精确罚函数的一个光滑近似[J]. 系统科学与数学, 2000, 20(4):499-505.

[7] 薛定宇,陈阳泉.高等应用数学问题的Matlab求解[M].北京:清华大学出版社,2004,37-42.

附录

1.障碍函数

function f=fun(x,r)

f=x(1,1)^2+x(2,1)^2-r*log(x(1,1)-1);

2.步长的函数

function f=fh(x0,h,s,r)

%h为步长

%s为方向

%r为1/t

x1=x0+h*s;

f=fun(x1,r);

3. 步长寻优函数

function h=fsearchh(x0,r,s)

%利用进退法确定高低高区间，利用黄金分割法进行求解

h1=0;%步长的初始点

st=0.001; %步长的步长

h2=h1+st;

f1=fh(x0,h1,s,r);

f2=fh(x0,h2,s,r);

if f1>f2

h3=h2+st;

f3=fh(x0,h3,s,r);

while f2>f3

h1=h2;

h2=h3;

h3=h3+st;

f2=f3;

f3=fh(x0,h3,s,r);

end

else

st=-st;

v=h1;

h1=h2;

h2=v;

v=f1;

f1=f2;

f2=v;

h3=h2+st;

f3=fh(x0,h3,s,r);

while f2>f3

h1=h2;

h2=h3;

h3=h3+st;

f2=f3;

f3=fh(x0,h3,s,r);

end

end

%得到高低高的区间

a=min(h1,h3);

b=max(h1,h3);

%利用黄金分割点法进行求解

h1=1+0.382*(b-a);

h2=1+0.618*(b-a);

f1=fh(x0,h1,s,r);

f2=fh(x0,h2,s,r);

while abs(a-b)>0.0001

if f1>f2

a=h1;

h1=h2;

f1=f2;

h2=a+0.618*(b-a);

f2=fh(x0,h2,s,r);

else

b=h2;

h2=h1;

f2=f1;

h1=a+0.382*(b-a);

f1=fh(x0,h1,s,r);

end

end

h=0.5*(a+b);

4. 迭代点的寻优函数

function f=fsearchx(x0,r,epson)

x00=x0;

m=length(x0);

s=zeros(m,1);

for i=1:m

s(i)=1;

h=fsearchh(x0,r,s);

x1=x0+h*s;

s(i)=0;

x0=x1;

end

while norm(x1-x00)>epson

x00=x1;

for i=1:m

s(i)=1;

h=fsearchh(x0,r,s);

x1=x0+h*s;

s(i)=0;

x0=x1;

end

end

f=x1;

5. 主程序

clear

clc

x0=[2;2]; %给定初始点

r=1;

c=0.1;

epson=0.001;

x1=fsearchx(x0,0.1,epson);

while norm(x0-x1)>epson

x0=x1;

r=r*c;

x1=fsearchx(x0,r,epson) ;

end

disp '函数的最优解为'

x1

6.算例函数图像

clc;clear;

xx=[1.01:0.01:2];

yy=[-2:0.01:2];

r=0.001;

[x,y]=meshgrid(xx,yy);

z1=x.^2+y.^2-r.*log2(x-1);

surf(x,y,z1),shading flat