matlab和python非線性規劃
不受約束的問題
非線性規劃問題分為兩種問題:無約束和約束。 不受約束的優化問題是形式上的問題
min
f
(
x
)
,
x
=
[
x
1
,
x
2
,
⋯
,
x
n
]
T
∈
R
n
(
1
)
\min \boldsymbol{f}\left( \boldsymbol{x} \right) ,\boldsymbol{x}=\left[ \boldsymbol{x}_1,\boldsymbol{x}_2,\cdots ,\boldsymbol{x}_{\boldsymbol{n}} \right] ^{\boldsymbol{T}}\in \mathbb{R}^{\boldsymbol{n}}\,\, \left( 1 \right)
一個一般的約束優化問題包括方程(1)中具有相等約束的目標函式:
E ( x ) = 0 ( 2 ) \boldsymbol{E}\left( \boldsymbol{x} \right) =0 \left( 2 \right) E(x)=0(2)
和不平等約束:
I ( x ) ⩽ 0 ( 3 ) \boldsymbol{I}\left( \boldsymbol{x} \right) \leqslant 0 \left( 3 \right) I(x)⩽0(3)
其中,
x
∈
R
n
,
E
(
x
)
∈
R
p
和
I
(
x
)
∈
R
q
\boldsymbol{x}\in \mathbb{R}^{\boldsymbol{n}},\boldsymbol{E}\left( \boldsymbol{x} \right) \in \mathbb{R}^{\boldsymbol{p}}和\,\,\boldsymbol{I}\left( \boldsymbol{x} \right) \in \mathbb{R}^{\boldsymbol{q}}
在本節中,將考慮上述問題(1)。如果 x ∗ = [ x 1 ∗ , x 2 ∗ , ⋯ , x n ∗ ] T ∈ R n \boldsymbol{x}^*=\left[ \boldsymbol{x}_{1}^{*},\boldsymbol{x}_{2}^{*},\cdots ,\boldsymbol{x}_{\boldsymbol{n}}^{*} \right] ^{\boldsymbol{T}}\in \mathbb{R}^{\boldsymbol{n}} x∗=[x1∗,x2∗,⋯,xn∗]T∈Rn是問題的最優解,則
g
(
x
∗
)
=
0
(
4
)
y
T
H
(
x
∗
)
y
⩾
0
,
∀
y
∈
R
n
(
5
)
\boldsymbol{g}\left( \boldsymbol{x}^* \right) =0 \left( 4 \right) \\\boldsymbol{y}^{\boldsymbol{T}}\boldsymbol{H}\left( \boldsymbol{x}^* \right) \boldsymbol{y}\geqslant 0,\forall \boldsymbol{y}\in \mathbb{R}^{\boldsymbol{n}}\,\, \left( 5 \right)
其中 g ( x ) = ∇ f ( x ) \boldsymbol{g}\left( \boldsymbol{x} \right) =\nabla \boldsymbol{f}\left( \boldsymbol{x} \right) g(x)=∇f(x) 是梯度向量, H ( x ) \boldsymbol{H}\left( \boldsymbol{x} \right) H(x) 是由下式定義的黑森矩陣,
g ( x ) = ∇ f ( x ) = [ ∂ f ( x ) ∂ x 1 ∂ f ( x ) x 2 ⋮ ∂ f ( x ) ∂ x n ] H ( x ) = [ ∂ 2 f ( x ) ∂ x 1 2 ∂ 2 f ( x ) ∂ x 1 ∂ x 2 ⋯ ∂ 2 f ( x ) ∂ x 1 ∂ x n ∂ 2 f ( x ) ∂ x 2 ∂ x 1 ∂ 2 f ( x ) ∂ x 2 2 ⋯ ∂ 2 f ( x ) ∂ x 2 ∂ x n ⋮ ⋮ ⋱ ⋮ ∂ 2 f ( x ) ∂ x n ∂ x 1 ∂ 2 f ( x ) ∂ x n ∂ x 2 ⋯ ∂ 2 f ( x ) ∂ x n 2 ] \boldsymbol{g}\left( \boldsymbol{x} \right) =\nabla \boldsymbol{f}\left( \boldsymbol{x} \right) =\left[ \begin{array}{c} \frac{\partial \boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_1}\\ \frac{\partial \boldsymbol{f}\left( \boldsymbol{x} \right)}{\boldsymbol{x}_2}\\ \vdots\\ \frac{\partial \boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{\boldsymbol{n}}}\\\end{array} \right] \\\boldsymbol{H}\left( \boldsymbol{x} \right) =\left[ \begin{matrix} \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{1}^{2}}& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_1\partial \boldsymbol{x}_2}& \cdots& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_1\partial \boldsymbol{x}_{\boldsymbol{n}}}\\ \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_2\partial \boldsymbol{x}_1}& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{2}^{2}}& \cdots& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_2\partial \boldsymbol{x}_{\boldsymbol{n}}}\\ \vdots& \vdots& \ddots& \vdots\\ \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{\boldsymbol{n}}\partial \boldsymbol{x}_1}& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{\boldsymbol{n}}\partial \boldsymbol{x}_2}& \cdots& \frac{\partial ^2\boldsymbol{f}\left( \boldsymbol{x} \right)}{\partial \boldsymbol{x}_{\boldsymbol{n}}^{2}}\\\end{matrix} \right] g(x)=∇f(x)=⎣⎢⎢⎢⎢⎡∂x1∂f(x)x2∂f(x)⋮∂xn∂f(x)⎦⎥⎥⎥⎥⎤H(x)=⎣⎢⎢⎢⎢⎢⎡∂x12∂2f(x)∂x2∂x1∂2f(x)⋮∂xn∂x1∂2f(x)∂x1∂x2∂2f(x)∂x22∂2f(x)⋮∂xn∂x2∂2f(x)⋯⋯⋱⋯∂x1∂xn∂2f(x)∂x2∂xn∂2f(x)⋮∂xn2∂2f(x)⎦⎥⎥⎥⎥⎥⎤
等式(5)表示黑森矩陣是一個正半定數。 式(4)-(5)是 x ∗ \boldsymbol{x}^* x∗的必要最優條件。 如果條件(5)替換為條件:
y T H ( x ∗ ) y > 0 , ∀ y ∈ R n ( 6 ) \boldsymbol{y}^{\boldsymbol{T}}\boldsymbol{H}\left( \boldsymbol{x}* \right) \boldsymbol{y}>0,\forall \boldsymbol{y}\in \mathbb{R}^{\boldsymbol{n}}\,\, \left( 6 \right) yTH(x∗)y>0,∀y∈Rn(6)
(即,黑森矩陣是正定的),則條件(4)-(6)是充分的最優性條件。
用於求解(1)的數值方法是迭代的。 他們將從解 x ∗ \boldsymbol{x}^* x∗的初始猜測 x ( 0 ) \boldsymbol{x}^{\left( 0 \right)} x(0)開始,然後構造解 x ( 0 ) , x ( 1 ) , x ( 2 ) , ⋯ , \boldsymbol{x}^{\left( 0 \right)},\boldsymbol{x}^{\left( 1 \right)},\boldsymbol{x}^{\left( 2 \right)},\cdots , x(0),x(1),x(2),⋯,的序列,其中 f ( x ( 0 ) ) > f ( x ( 1 ) ) > ⋯ \boldsymbol{f}\left( \boldsymbol{x}^{\left( 0 \right)} \right) >\boldsymbol{f}\left( \boldsymbol{x}^{\left( 1 \right)} \right) >\cdots f(x(0))>f(x(1))>⋯ 在迭代 k \boldsymbol{k} k時,數值優化中從點 x ( k ) \boldsymbol{x}^{\left( \boldsymbol{k} \right)} x(k)移至點 x ( k + 1 ) \boldsymbol{x}^{\left( \boldsymbol{k}+1 \right)} x(k+1)通常需要經歷兩個步驟:第一步是搜尋方向 p ( k ) \boldsymbol{p}^{\left( \boldsymbol{k} \right)} p(k)確定,然後在第二步中沿 p ( k ) \boldsymbol{p}^{\left( \boldsymbol{k} \right)} p(k)方向進行線搜尋(其中 ∥ p ( k + 1 ) ∥ = 1 \left\| \boldsymbol{p}^{\left( \boldsymbol{k}+1 \right)} \right\| =1 ∥∥p(k+1)∥∥=1)以找到最小點 x ( k + 1 ) = x ( k ) + α p ( k ) \boldsymbol{x}^{\left( \boldsymbol{k}+1 \right)}=\boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha p}^{\left( \boldsymbol{k} \right)} x(k+1)=x(k)+αp(k),使得 f ( x ( k + 1 ) ) < f ( x ( k ) ) \boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k}+1 \right)} \right) <\boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)} \right) f(x(k+1))<f(x(k)) 。 如果滿足以下條件,則迭代過程在解 x ( k ) \boldsymbol{x}^{\left( \boldsymbol{k} \right)} x(k)處停止。
∥ g ( x ( k ) ) ∥ < ε \left\| \boldsymbol{g}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)} \right) \right\| <\boldsymbol{\varepsilon } ∥∥∥g(x(k))∥∥∥<ε
其中 ε < 0 \boldsymbol{\varepsilon }<0 ε<0是任意小的正實數,而 H ( x ( k ) ) \boldsymbol{H}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)} \right) H(x(k)) 是正半定數。
無約束優化的數值優化技術關注兩個問題:
-
線搜尋問題:給定一個函式 f ( x ) \boldsymbol{f}\left( \boldsymbol{x} \right) f(x) ,其梯度 g ( x ) = ∇ f ( x ) \boldsymbol{g}\left( \boldsymbol{x} \right) =\nabla \boldsymbol{f}\left( \boldsymbol{x} \right) g(x)=∇f(x) ,在迭代 k \boldsymbol{k} k下降方向 p ( k ) \boldsymbol{p}^{\left( \boldsymbol{k} \right)} p(k)上的一個點 x ( k ) \boldsymbol{x}^{\left( \boldsymbol{k} \right)} x(k),找到 α ( k ) > 0 \boldsymbol{\alpha }^{\left( \boldsymbol{k} \right)}>0 α(k)>0使得對於 α = α ( k ) \boldsymbol{\alpha }=\boldsymbol{\alpha }^{\left( \boldsymbol{k} \right)} α=α(k),沿著射線 x ( k ) + α p ( k ) \boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha p}^{\left( \boldsymbol{k} \right)} x(k)+αp(k)函式 f \boldsymbol{f} f最小化,即
α ( k ) = a r g min α > 0 { f ( x ( k ) + α p ( k ) ) } \boldsymbol{\alpha }^{\left( \boldsymbol{k} \right)}=\boldsymbol{arg}\underset{\boldsymbol{\alpha }>0}{\min}\left\{ \boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha p}^{\left( \boldsymbol{k} \right)} \right) \right\} α(k)=argα>0min{f(x(k)+αp(k))}
-
搜尋方向問題:給定一個函式 f ( x ) \boldsymbol{f}\left( \boldsymbol{x} \right) f(x) 和點 x ( k ) \boldsymbol{x}^{\left( \boldsymbol{k} \right)} x(k),在迭代 k \boldsymbol{k} k處,找到一個單位向量 p ( k ) \boldsymbol{p}^{\left( \boldsymbol{k} \right)} p(k),使得 f ( x ( k ) + α p ( k ) ) \boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha p}^{\left( \boldsymbol{k} \right)} \right) f(x(k)+αp(k)) 對於 0 < α < α max 0<\boldsymbol{\alpha }<\boldsymbol{\alpha }_{\max} 0<α<αmax是遞減函式。 即 p ( k ) \boldsymbol{p}^{\left( \boldsymbol{k} \right)} p(k)是下降方向。
數值優化方法的不同之處在於確定搜尋方向以及對梯度向量和黑森矩陣進行近似和(或)更新的方法。
線搜尋演算法
線搜尋演算法有兩種:精確線搜尋和不精確線搜尋演算法。 精確線搜尋演算法會尋找精確步長 α ( α = α E x a c t ) \boldsymbol{\alpha }\left( \boldsymbol{\alpha }=\boldsymbol{\alpha }_{\boldsymbol{Exact}} \right) α(α=αExact) ,以便
α E x a c t = a r g min α > 0 { f ( x ( k ) + α p ( k ) ) } \boldsymbol{\alpha }_{\boldsymbol{Exact}}=\boldsymbol{arg}\underset{\boldsymbol{\alpha }>0}{\min}\left\{ \boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha p}^{\left( \boldsymbol{k} \right)} \right) \right\} αExact=argα>0min{f(x(k)+αp(k))}
在數字上不切實際。 不精確的線搜尋演算法尋找與精確步長 α E x a c t \boldsymbol{\alpha }_{\boldsymbol{Exact}} αExact近似的值。 從初始猜測 α 0 \boldsymbol{\alpha }^0 α0開始,逐步構建步長序列 α 0 , α 1 , α 2 , ⋯ , α k \boldsymbol{\alpha }^0,\boldsymbol{\alpha }^1,\boldsymbol{\alpha }^2,\cdots ,\boldsymbol{\alpha }^{\boldsymbol{k}} α0,α1,α2,⋯,αk的序列,以使 α k ≈ α E x a c t \boldsymbol{\alpha }^{\boldsymbol{k}}\approx \boldsymbol{\alpha }_{\boldsymbol{Exact}} αk≈αExact迭代地完成。
最著名的不精確線搜尋演算法是回溯線搜尋演算法,該演算法使用兩個引數 a \boldsymbol{a} a和 b \boldsymbol{b} b並基於Armijo條件:
f ( x ( k ) + α k p ( k ) ) < f ( x ) + a α k g ( x ( k ) ) T p ( k ) \boldsymbol{f}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)}+\boldsymbol{\alpha }^{\boldsymbol{k}}\boldsymbol{p}^{\left( \boldsymbol{k} \right)} \right) <\boldsymbol{f}\left( \boldsymbol{x} \right) +\boldsymbol{a\alpha }^{\boldsymbol{k}}\boldsymbol{g}\left( \boldsymbol{x}^{\left( \boldsymbol{k} \right)} \right) ^{\boldsymbol{T}}\boldsymbol{p}^{\left( \boldsymbol{k} \right)} f(x(k)+αkp(k))<f(x)+aαkg(x(k))Tp(k)
引數 a \boldsymbol{a} a與Armijo的終止條件相關。 引數 0 < b < 1 0<\boldsymbol{b}<1 0<b<1用於在每次迭代中更新步長,其中 α j = b α j − 1 \boldsymbol{\alpha }^{\boldsymbol{j}}=\boldsymbol{b\alpha }^{\boldsymbol{j}-1} αj=bαj−1,從大步長 α 0 \boldsymbol{\alpha }^0 α0(通常= 1)開始。
MATLAB函式LineSearch.m接收一個函式 f \boldsymbol{f} f,代表起始向量 x \boldsymbol{x} x的 f \boldsymbol{f} f梯度的向量 g \boldsymbol{g} g和單位(方向)向量 p \boldsymbol{p} p並返回最佳步長 α \boldsymbol{\alpha } α:
function alpha = LineSearch(f, g, x, p)
a = 0.3 ; b = 0.9 ;
alpha = 1.0 ;
while f(x+alpha*p) > f(x) + a*alpha*g(x)'*p
alpha = b*alpha ;
end
end
函式LineSearch.py的Python程式碼是:
import numpy as np
def LineSearch(f, g, x, p):
a, b = 0.3, 0.9
alpha = 1.0
while f(x+alpha*p) > f(x) + a*alpha*np.dot(g(x), p):
alpha *= b
return alpha