1. 程式人生 > 其它 >matlab和python非線性規劃

matlab和python非線性規劃

技術標籤:Pythonmatlabpython

Python MATLAB

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)

minf(x),x=[x1,x2,,xn]TRn(1)

一個一般的約束優化問題包括方程(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}}

xRn,E(x)RpI(x)Rq

在本節中,將考慮上述問題(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]TRn是問題的最優解,則

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)=0(4)yTH(x)y0,yRn(5)

其中 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)=x1f(x)x2f(x)xnf(x)H(x)=x122f(x)x2x12f(x)xnx12f(x)x1x22f(x)x222f(x)xnx22f(x)x1xn2f(x)x2xn2f(x)xn22f(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,yRn(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)) 是正半定數。

無約束優化的數值優化技術關注兩個問題:

  1. 線搜尋問題:給定一個函式 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))}

  2. 搜尋方向問題:給定一個函式 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αj1,從大步長 α 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

最陡下降法

牛頓法

擬牛頓法

用MATLAB解決無約束的優化問題

用Gekko解決無約束的優化問題

解決約束的優化問題

MATLAB fmincon函式解決約束優化問題

Python解決約束最小化問題

Gekko Python解決約束優化

詳情參閱http://viadean.com/matlab_python_nonlinear_prg.html