最优化方法Python计算:求解约束优化问题的拉格朗日乘子算法

从仅有等式约束的问题入手。设优化问题(7.8)
{ minimize f ( x ) s.t. h ( x ) = o ( 1 ) \begin{cases} \text{minimize}\quad\quad f(\boldsymbol{x})\\ \text{s.t.}\quad\quad\quad \boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o} \end{cases}\quad\quad(1) {minimizef(x)s.t.h(x)=o(1)
f ( x ) f(\boldsymbol{x}) f(x) h ( x ) \boldsymbol{h}(\boldsymbol{x}) h(x)均在 R n \text{R}^n Rn上二阶连续可微,且 ( x 0 λ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix} (x0λ0)是问题(1)满足二阶充分条件的KKT点。即 ∇ f ( x 0 ) − ∇ h ( x 0 ) λ 0 = 0 \nabla f(\boldsymbol{x}_0)-\nabla\boldsymbol{h}(\boldsymbol{x}_0)\boldsymbol{\lambda}_0=0 f(x0)h(x0)λ0=0 ∀ d ∈ { d ∣ d ≠ o , ∇ h ( x ) ⊤ d = o } \forall\boldsymbol{d}\in\{\boldsymbol{d}|\boldsymbol{d}\not=\boldsymbol{o},\nabla\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{d}=\boldsymbol{o}\} d{dd=o,h(x)d=o} d ⊤ ∇ x x L ( x 0 , λ 0 ) d > 0 \boldsymbol{d}^\top\nabla_{\boldsymbol{xx}} L(\boldsymbol{x}_0,\boldsymbol{\lambda}_0)\boldsymbol{d}>0 dxxL(x0,λ0)d>0。其中, L ( x , λ ) = f ( x ) − λ ⊤ h ( x ) L(\boldsymbol{x},\boldsymbol{\lambda})=f(\boldsymbol{x})-\boldsymbol{\lambda}^\top\boldsymbol{h}(\boldsymbol{x}) L(x,λ)=f(x)λh(x)为问题(1)的拉格朗日函数,其中的 λ ∈ R l \boldsymbol{\lambda}\in\text{R}^l λRl为拉格朗日乘子。
定义问题(7.8)的{\heiti{增广拉格朗日函数}}
L ( x , λ , σ ) = f ( x ) − λ ⊤ h ( x ) + σ 2 h ( x ) ⊤ h ( x ) ( 2 ) L(\boldsymbol{x},\boldsymbol{\lambda},\sigma)=f(\boldsymbol{x})-\boldsymbol{\lambda}^\top\boldsymbol{h}(\boldsymbol{x})+\frac{\sigma}{2}\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{h}(\boldsymbol{x})\quad\quad{(2)} L(x,λ,σ)=f(x)λh(x)+2σh(x)h(x)(2)
其中,罚因子 σ > 0 \sigma>0 σ>0。与问题(1)的拉格朗日函数 L ( x , λ ) L(\boldsymbol{x},\boldsymbol{\lambda}) L(x,λ)相比,增广拉格朗日函数 L ( x , λ , σ ) L(\boldsymbol{x},\boldsymbol{\lambda},\sigma) L(x,λ,σ)多了个罚项 σ 2 h ( x ) ⊤ h ( x ) \frac{\sigma}{2}\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{h}(\boldsymbol{x}) 2σh(x)h(x)。对问题(1)的增广拉格朗日函数(2),不难验证 ∀ x ∈ Ω \forall\boldsymbol{x}\in\Omega xΩ L ( x , λ , σ ) = f ( x ) L(\boldsymbol{x},\boldsymbol{\lambda},\sigma)=f(\boldsymbol{x}) L(x,λ,σ)=f(x)。可以证明以下命题:
定理1 ( x 0 λ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix} (x0λ0)是问题(1)满足二阶充分条件的KKT点,则存在 σ ′ ≥ 0 \sigma'\geq0 σ0,使得对任意的 σ > σ ′ \sigma>\sigma' σ>σ x 0 = arg ⁡ min ⁡ x ∈ R n L ( x , λ 0 , σ ) . \boldsymbol{x}_0=\arg\min_{\boldsymbol{x}\in\text{R}^n}L(\boldsymbol{x},\boldsymbol{\lambda}_0,\sigma). x0=argxRnminL(x,λ0,σ).
反之,若 x 0 = arg ⁡ min ⁡ x ∈ R n L ( x , λ 0 , σ ) \boldsymbol{x}_0=\arg\min\limits_{\boldsymbol{x}\in\text{R}^n}L(\boldsymbol{x},\boldsymbol{\lambda}_0,\sigma) x0=argxRnminL(x,λ0,σ),且 h ( x 0 ) = o \boldsymbol{h}(\boldsymbol{x}_0)=\boldsymbol{o} h(x0)=o,则 x 0 \boldsymbol{x}_0 x0是问题(1)的局部最优解,即 x 0 = arg ⁡ min ⁡ x ∈ Ω f ( x ) \boldsymbol{x}_0=\arg\min\limits_{\boldsymbol{x}\in\Omega}f(\boldsymbol{x}) x0=argxΩminf(x)
定理1将存在最优解满足二阶充分条件KKT点 ( x 0 λ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix} (x0λ0)的等式约束优化问题(1)求解,转化为对其增广拉格朗日函数 L ( x , λ 0 , σ ) L(\boldsymbol{x},\boldsymbol{\lambda}_0,\sigma) L(x,λ0,σ)作为目标函数的辅助无约束优化问题的求解,其中 σ > 0 \sigma>0 σ>0是一个充分大的实数。然而,一般情况下,拉格朗日乘子 λ 0 \boldsymbol{\lambda}_0 λ0未必已知。需设法用一系列 { λ k } \{\boldsymbol{\lambda}_k\} {λk}去逼近 λ 0 \boldsymbol{\lambda}_0 λ0
利用定理1,我们得到由初始的乘子 λ 1 \boldsymbol{\lambda}_1 λ1和罚因子 σ 1 > 0 \sigma_1>0 σ1>0出发,计算问题(1)的满足二阶充分条件的KKT点 ( x 0 λ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix} (x0λ0)近似值的迭代方法:计算
{ x k + 1 = arg ⁡ min ⁡ x ∈ R n L ( x , λ k , σ k ) λ k + 1 = λ k − σ k h ( x k ) . \begin{cases} \boldsymbol{x}_{k+1}=\arg\min\limits_{\boldsymbol{x}\in\text{R}^n}L(\boldsymbol{x},\boldsymbol{\lambda}_k,\sigma_k)\\ \boldsymbol{\lambda}_{k+1}=\boldsymbol{\lambda}_k-\sigma_k\boldsymbol{h}(\boldsymbol{x}_{k}) \end{cases}. {xk+1=argxRnminL(x,λk,σk)λk+1=λkσkh(xk).
对给定的容错误差 ε > 0 \varepsilon>0 ε>0,若 ∥ h ( x k + 1 ) ∥ < ε \lVert\boldsymbol{h}(\boldsymbol{x}_{k+1})\rVert<\varepsilon h(xk+1)∥<ε,则停止迭代。 ( x k + 1 λ k + 1 ) \begin{pmatrix}\boldsymbol{x}_{k+1}\\\boldsymbol{\lambda}_{k+1}\end{pmatrix} (xk+1λk+1)即为 ( x 0 λ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix} (x0λ0)近似值。\par
否则,即 ∥ h ( x k + 1 ) ∥ ≥ ε \lVert\boldsymbol{h}(\boldsymbol{x}_{k+1})\rVert\geq\varepsilon h(xk+1)∥ε,用给定的放大系数 c > 1 c>1 c>1计算
σ k + 1 = { σ k ∥ h k + 1 ∥ ∥ h k ∥ < 1 c × σ k ∥ h k + 1 ∥ ∥ h k ∥ ≥ 1 , \sigma_{k+1}=\begin{cases} \sigma_k&\frac{\lVert\boldsymbol{h}_{k+1}\rVert}{\lVert\boldsymbol{h}_{k}\rVert}<1\\ c\times\sigma_k&\frac{\lVert\boldsymbol{h}_{k+1}\rVert}{\lVert\boldsymbol{h}_{k}\rVert}\geq1 \end{cases}, σk+1={σkc×σkhkhk+1<1hkhk+11,
进入下一轮迭代。
以上讨论的对等式约束优化问题的乘子算法是由 Powell和Hestenes于1969年提出的,常称为PH算法
对仅有不等式约束
{ minimize f ( x ) s.t.   g ( x ) ≥ o \begin{cases} \text{minimize}\quad f(\boldsymbol{x})\\ \text{s.t.\ \ }\quad\quad\boldsymbol{g}(\boldsymbol{x})\geq\boldsymbol{o} \end{cases} {minimizef(x)s.t.  g(x)o
其中, f : R n → R f:\text{R}^n\rightarrow\text{R} f:RnR g : R n → R m \boldsymbol{g}:\text{R}^n\rightarrow\text{R}^m g:RnRm R n \text{R}^n Rn上二阶连续可微。可行域 Ω = { x ∣ g ( x ) ≥ o } \Omega=\{\boldsymbol{x}|\boldsymbol{g}(\boldsymbol{x})\geq\boldsymbol{o}\} Ω={xg(x)o}。添加松弛变量 s ≥ o \boldsymbol{s}\geq\boldsymbol{o} so,构造等价问题
{ minimize f ( x ) s.t.   g ( x ) − s = o s ≥ o \begin{cases} \text{minimize}\quad f(\boldsymbol{x})\\ \text{s.t.\ \ }\quad\quad\boldsymbol{g}(\boldsymbol{x})-\boldsymbol{s}=\boldsymbol{o}\\ \quad\quad\quad\quad\boldsymbol{s}\geq\boldsymbol{o} \end{cases} minimizef(x)s.t.  g(x)s=oso
通过变换,上述问题可转换成等价的
{ minimize f ( x ) s.t. h 1 ( x ) = o . \begin{cases} \text{minimize}\quad\quad f(\boldsymbol{x})\\ \text{s.t.}\quad\quad\quad \boldsymbol{h}_1(\boldsymbol{x})=\boldsymbol{o} \end{cases}. {minimizef(x)s.t.h1(x)=o.
其中, h 1 ( x ) = min ⁡ ( 1 σ μ , g ( x ) ) \boldsymbol{h}_1(\boldsymbol{x})=\boldsymbol{\min}\left(\frac{1}{\sigma}\boldsymbol{\mu},\boldsymbol{g}(\boldsymbol{x})\right) h1(x)=min(σ1μ,g(x))
相仿地,对一般的优化问题
{ minimize f ( x ) s.t.   h ( x ) = o g ( x ) ≥ o ( 3 ) \begin{cases} \text{minimize}\quad f(\boldsymbol{x})\\ \text{s.t.\ \ }\quad\quad\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o}\\ \quad\quad\quad\quad\boldsymbol{g}(\boldsymbol{x})\geq\boldsymbol{o} \end{cases}\quad\quad(3) minimizef(x)s.t.  h(x)=og(x)o(3)
其中 f ( x ) f(\boldsymbol{x}) f(x) h ( x ) \boldsymbol{h}(\boldsymbol{x}) h(x) g ( x ) \boldsymbol{g}(\boldsymbol{x}) g(x) R n \text{R}^n Rn上二阶连续可微,可行域 Ω = { x ∣ h ( x ) = o , g ( x ≥ o ) } \Omega=\{\boldsymbol{x}|\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o},\boldsymbol{g}(\boldsymbol{x}\geq\boldsymbol{o})\} Ω={xh(x)=o,g(xo)}。为求解其满足二阶充分条件的KKT点 ( x 0 λ 0 μ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\\\boldsymbol{\mu}_0\end{pmatrix} x0λ0μ0 ,将其转换为与之等价的问题
{ minimize f ( x ) s.t.   h ( x ) = o h 1 ( x ) = o . ( 4 ) \begin{cases} \text{minimize}\quad f(\boldsymbol{x})\\ \text{s.t.\ \ }\quad\quad\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o}\\ \quad\quad\quad\quad\boldsymbol{h}_1(\boldsymbol{x})=\boldsymbol{o} \end{cases}.\quad\quad{(4)} minimizef(x)s.t.  h(x)=oh1(x)=o.(4)
其增广拉格朗日函数为
L ( x , λ , μ , σ ) = f ( x ) − λ ⊤ h ( x ) − μ ⊤ h 1 ( x ) + σ 2 ( h ( x ) ⊤ h ( x ) + h 1 ( x ) ⊤ h 1 ( x ) ) L(\boldsymbol{x},\boldsymbol{\lambda},\boldsymbol{\mu},\sigma)=f(\boldsymbol{x})-\boldsymbol{\lambda}^\top\boldsymbol{h}(x)-\boldsymbol{\mu}^\top\boldsymbol{h_1}(\boldsymbol{x})+\frac{\sigma}{2}(\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{h}(\boldsymbol{x})+\boldsymbol{h}_1(\boldsymbol{x})^\top\boldsymbol{h}_1(\boldsymbol{x})) L(x,λ,μ,σ)=f(x)λh(x)μh1(x)+2σ(h(x)h(x)+h1(x)h1(x))
对初始乘子 λ 1 \boldsymbol{\lambda}_1 λ1 μ 1 \boldsymbol{\mu}_1 μ1,罚因子 σ 1 \sigma_1 σ1,用拉格朗日乘子法计算问题(7.10)满足二阶充分条件的KKT点 ( x 0 λ 0 μ 0 ) \begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\\\boldsymbol{\mu}_0\end{pmatrix} x0λ0μ0 的近似值,迭代式为
{ x k + 1 = arg ⁡ min ⁡ x ∈ R n L ( x , λ k , μ k , σ k ) λ k + 1 = λ k − σ k h ( x k ) μ k + 1 = max ⁡ ( o , μ k − σ k g ( x k ) ) , \begin{cases} \boldsymbol{x}_{k+1}=\arg\min\limits_{\boldsymbol{x}\in\text{R}^n}L(\boldsymbol{x},\boldsymbol{\lambda}_k,\boldsymbol{\mu}_k,\sigma_k)\\ \boldsymbol{\lambda}_{k+1}=\boldsymbol{\lambda}_k-\sigma_k\boldsymbol{h}(\boldsymbol{x}_k)\\ \boldsymbol{\mu}_{k+1}=\boldsymbol{\max}(\boldsymbol{o},\boldsymbol{\mu}_k-\sigma_k\boldsymbol{g}(\boldsymbol{x}_k)) \end{cases}, xk+1=argxRnminL(x,λk,μk,σk)λk+1=λkσkh(xk)μk+1=max(o,μkσkg(xk)),
对给定的容错误差 ε > 0 \varepsilon>0 ε>0,记 β = ∥ h ( x k ) ∥ + ∥ h 1 ( x k ) ∥ \beta=\lVert\boldsymbol{h}(\boldsymbol{x}_k)\rVert+\lVert\boldsymbol{h}_1(\boldsymbol{x}_k)\rVert β=h(xk)∥+h1(xk)∥,则迭代终止条件为
β < ε , \beta<\varepsilon, β<ε,
罚因子扩大条件为
β n e w β o l d = ∥ h ( x k + 1 ) ∥ + ∥ h 1 ( x k + 1 ) ∥ ∥ h ( x k ) ∥ + ∥ h 1 ( x k ) ∥ ≥ η . \frac{\beta_{new}}{\beta_{old}}=\frac{\lVert\boldsymbol{h}(\boldsymbol{x}_{k+1})\rVert+\lVert\boldsymbol{h}_1(\boldsymbol{x}_{k+1})\rVert}{\lVert\boldsymbol{h}(\boldsymbol{x}_k)\rVert+\lVert\boldsymbol{h}_1(\boldsymbol{x}_k)\rVert}\geq\eta. βoldβnew=h(xk)∥+h1(xk)∥h(xk+1)∥+h1(xk+1)∥η.
其中 η ∈ ( 0 , 1 ) \eta\in(0,1) η(0,1)。这一算法是Rockfellar于1973年,在PH算法的基础上提出的,常称为PHR算法。下列代码实现PHR算法。

import numpy as np											#导入numpy
from scipy.optimize import minimize, OptimizeResult			#导入minimize和OptimizeResult
def phr(f, x1, h = None, g = None, eps = 1e-5):k = 0sigmak = 10c = 2.5eta = 0.8if callable(h):											#有等式约束l = h(x1).sizelamdak = 0.1 * np.ones(l)else:													#无等式约束lamdak = np.zeros(1)h = lambda x: np.zeros(1)if callable(g):											#有不等式约束m = g(x1).sizemuk = 0.1 * np.ones(m)h1 = lambda x: np.minimum(muk / sigmak, g(x))else:													#无不等式约束muk = np.zeros(1)h1 = g = lambda x: np.zeros(1)m = 1zero = np.zeros(m)										#零向量P = lambda x: np.dot(h(x), h(x)) + np.dot(h1(x), h1(x))	#罚函数L = lambda x: f(x) - np.dot(lamdak, h(x)) - \			#增广拉格朗日函数np.dot(muk, h1(x)) + 0.5 * sigmak * P(x)xk = x1													#迭代点bnew = bold = np.linalg.norm(h(xk)) + np.linalg.norm(h1(xk))#新、旧beta值while bnew >= eps:k += 1xk = minimize(L, xk).x								#更新迭代点lamdak = lamdak - sigmak * h(xk)					#更新乘子lambdamuk = np.maximum(zero, muk - sigmak * g(xk))		#更新乘子muif bnew / bold >= eta:								#须扩大罚因子sigmak *= cbold = bnew											#保存betabnew = np.linalg.norm(h(xk)) + np.linalg.norm(h1(xk))#更新betareturn OptimizeResult(fun = f(xk), x = xk, nit = k)

程序的第3~37行定义的函数phr实现图7.4所示的乘子算法PHR。该函数的5个参数f,x1,h,g和eps分别表示优化问题(3)的目标函数 f ( x ) f(\boldsymbol{x}) f(x),初始迭代点 x 1 \boldsymbol{x}_1 x1,等式约束函数 h ( x ) \boldsymbol{h}(\boldsymbol{x}) h(x),不等式约束函数 h ( x ) \boldsymbol{h}(\boldsymbol{x}) h(x)和容错误差 ε \varepsilon ε。其中,h,g和eps的缺省值分别为None、None和 1 0 − 5 10^{-5} 105
函数体内,第4~7行分别初始化迭代次数 k k k,罚因子 σ k \sigma_k σk,罚因子放大系数 c c c和罚因子放大判别值 η \eta η
第8~13行的if-else语句根据是否有等式约束,初始化乘子 λ \boldsymbol{\lambda} λ:若是,初始化为含有 l l l个0.1的向量,否则置为0。相仿地,第14~21行的if-else语句设置乘子 μ \boldsymbol{\mu} μ函数及 h 1 ( x ) = min ⁡ { μ σ k , g ( x ) } \boldsymbol{h}_1(\boldsymbol{x})=\min\left\{\frac{\boldsymbol{\mu}}{\sigma_k},\boldsymbol{g}(\boldsymbol{x})\right\} h1(x)=min{σkμ,g(x)}
第23~25行定义罚函数和增广拉格朗日函数
P ( x ) = h ( x ) ⊤ h ( x ) + h 1 ( x ) ⊤ h 1 ( x ) L ( x ) = f ( x ) − λ ⊤ h ( x ) − μ ⊤ h 1 ( x ) + σ k 2 P ( x ) . \begin{array}{l} P(\boldsymbol{x})=\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{h}(\boldsymbol{x})+\boldsymbol{h}_1(\boldsymbol{x})^\top\boldsymbol{h}_1(\boldsymbol{x})\\ L(\boldsymbol{x})=f(\boldsymbol{x})-\boldsymbol{\lambda}^\top\boldsymbol{h}(x)-\boldsymbol{\mu}^\top\boldsymbol{h}_1(\boldsymbol{x})+\frac{\sigma_k}{2}P(\boldsymbol{x}) \end{array}. P(x)=h(x)h(x)+h1(x)h1(x)L(x)=f(x)λh(x)μh1(x)+2σkP(x).
第26行初始化迭代点xk,第27行将 β n e w \beta_{new} βnew β o l d \beta_{old} βold初始化为 ∥ h ( x k ) ∥ + ∥ h 1 ( x k ) ∥ \lVert\boldsymbol{h}(\boldsymbol{x}_k)\rVert+\lVert\boldsymbol{h}_1(\boldsymbol{x}_k)\rVert h(xk)∥+h1(xk)∥
第28~36行的while循环执行算法的迭代。其中,第30~32行按
{ x k + 1 = arg ⁡ min ⁡ x ∈ R n L ( x , λ k , μ k , σ k ) λ k + 1 = λ k − σ k h ( x k ) μ k + 1 = max ⁡ ( o , μ k − σ k g ( x k ) ) , \begin{cases} \boldsymbol{x}_{k+1}=\arg\min\limits_{\boldsymbol{x}\in\text{R}^n}L(\boldsymbol{x},\boldsymbol{\lambda}_k,\boldsymbol{\mu}_k,\sigma_k)\\ \boldsymbol{\lambda}_{k+1}=\boldsymbol{\lambda}_k-\sigma_k\boldsymbol{h}(\boldsymbol{x}_k)\\ \boldsymbol{\mu}_{k+1}=\boldsymbol{\max}(\boldsymbol{o},\boldsymbol{\mu}_k-\sigma_k\boldsymbol{g}(\boldsymbol{x}_k)) \end{cases}, xk+1=argxRnminL(x,λk,μk,σk)λk+1=λkσkh(xk)μk+1=max(o,μkσkg(xk)),
更新迭代点xk,乘子lamdk和muk。第33~34行根据是否
β n e w β o l d ≥ η \frac{\beta_{new}}{\beta_{old}}\geq\eta βoldβnewη
决定扩大罚因子 σ k \sigma_k σk。第35、36行更新 β o l d \beta_{old} βold β n e w \beta_{new} βnew,以备进入下一轮迭代。一旦测得
β n e w < ε \beta_{new}<\varepsilon βnew<ε
迭代结束,第37行返回最优解近似值xk,最优值近似值f(xk)和迭代次数k。
例1:用程序7.11定义的phr函数求解优化问题
{ minimize f ( x ) = ( x 1 − 2 ) 2 + ( x 2 − 1 ) 2 s.t.   x 1 − 2 x 2 + 1 = 0 1 4 x 1 2 + x 2 2 ≤ 1 , \begin{cases} \text{minimize}\quad f(\boldsymbol{x})=(x_1-2)^2+(x_2-1)^2\\ \text{s.t.\ \ }\quad\quad\quad x_1-2x_2+1=0\\ \quad\quad\quad\quad\quad \frac{1}{4}x_1^2+x_2^2\leq1 \end{cases}, minimizef(x)=(x12)2+(x21)2s.t.  x12x2+1=041x12+x221,
给定初始迭代点 x 1 = ( 3 3 ) \boldsymbol{x}_1=\begin{pmatrix}3\\3\end{pmatrix} x1=(33)
:本题中目标函数 f ( x ) = ( x 1 − 2 ) 2 + ( x 2 − 1 ) 2 f(\boldsymbol{x})=(x_1-2)^2+(x_2-1)^2 f(x)=(x12)2+(x21)2,等式约束函数 h ( x ) = x 1 − 2 x 2 + 1 h(\boldsymbol{x})=x_1-2x_2+1 h(x)=x12x2+1,不等式约束函数 g ( x ) = 1 − 1 4 x 1 2 − x 2 2 g(\boldsymbol{x})=1-\frac{1}{4}x_1^2-x_2^2 g(x)=141x12x22。下列代码完成计算。

import numpy as np								#导入numpy
from utility import phr							#导入phr
np.set_printoptions(precision=4)				#设置输出精度
f = lambda x: (x[0] - 2) ** 2 + (x[1] - 1) ** 2	#目标函数
h = lambda x: x[0] - 2 * x[1] + 1				#等式约束函数
g = lambda x: 1 - 0.25 * x[0] ** 2 - x[1] ** 2	#不等式约束函数
x1 = np.array([3, 3])							#初始迭代点
print(phr(f, x1, h, g))

利用代码内注释信息不难理解程序。运行程序,输出

 fun: 1.3934640206427056nit: 5x: array([0.8229, 0.9114])

意为phr用5次迭代,以 1 0 − 5 10^{-5} 105的容错误差,算得最优解近似值 x 0 ≈ ( 0.8229 0.9114 ) \boldsymbol{x}_0\approx\begin{pmatrix}0.8229\\0.9114\end{pmatrix} x0(0.82290.9114),最优值的近似值 f ( x 0 ) ≈ 1.3934 f(\boldsymbol{x}_0)\approx1.3934 f(x0)1.3934
例2:用phr函数求解求解线性规划
{ minimize x 1 − x 2 s.t.      − x 1 + 2 x 2 + x 3 ≤ 2 − 4 x 1 + 4 x 2 − x 3 = 4 x 1 − x 3 = 0 x 1 , x 2 , x 3 ≥ 0 , \begin{cases} \text{minimize}\quad\quad x_1-x_2\\ \text{s.t.\ \ \ \ \ }\quad\quad -x_1+2x_2+x_3\leq2\\ \quad\quad\quad\quad\quad -4x_1+4x_2-x_3=4\\ \quad\quad\quad\quad\quad x_1-x_3=0\\ \quad\quad\quad\quad\quad x_1,x_2,x_3\geq0 \end{cases}, minimizex1x2s.t.     x1+2x2+x324x1+4x2x3=4x1x3=0x1,x2,x30,
给定初始迭代点 x 1 = o \boldsymbol{x}_1=\boldsymbol{o} x1=o。\par
:这个线性规划的数据矩阵为
c = ( 1 − 1 0 ) , A e q = ( − 4 4 − 1 1 0 − 1 ) , b e q = ( 4 0 ) , A i q = ( 1 − 2 − 1 1 0 0 0 1 0 0 0 1 ) , b i q = ( − 2 0 0 0 ) . \boldsymbol{c}=\begin{pmatrix}1\\-1\\0\end{pmatrix},\boldsymbol{A}_{eq}=\begin{pmatrix}-4&4&-1\\1&0&-1\end{pmatrix},\boldsymbol{b}_{eq}=\begin{pmatrix}4\\0\end{pmatrix},\boldsymbol{A}_{iq}=\begin{pmatrix}1&-2&-1\\1&0&0\\0&1&0\\0&0&1\end{pmatrix},\boldsymbol{b}_{iq}=\begin{pmatrix}-2\\0\\0\\0\end{pmatrix}. c= 110 ,Aeq=(414011),beq=(40),Aiq= 110020101001 ,biq= 2000 .
于是
f ( x ) = c ⊤ x , h ( x ) = A e q − b e q , g ( x ) = A i q − b i p . f(\boldsymbol{x})=\boldsymbol{c}^\top\boldsymbol{x},\quad\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{A}_{eq}-\boldsymbol{b}_{eq},\quad\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{A}_{iq}-\boldsymbol{b}_{ip}. f(x)=cx,h(x)=Aeqbeq,g(x)=Aiqbip.,下列程序完成计算。

import numpy as np					#导入numpy
c = np.array([1, -1, 0])			#向量c
Ae = np.array([[-4, 4, -1],			#矩阵Aeq[1, 0, -1]])
be = np.array([4, 0])				#向量beq
Ai = np.array([[1, -2, -1],			#矩阵Aiq[1, 0, 0],[0, 1, 0],[0, 0, 1]])
bi = np.array([-2, 0, 0, 0])		#向量biq
f = lambda x: np.dot(c, x)			#目标函数
h = lambda x: np.matmul(Ae, x) - be	#等式约束函数
g = lambda x: np.matmul(Ai, x) - bi	#不等式约束函数
x1 = np.zeros(3)					#初始迭代点
print(phr(f, x1, h, g))

运行程序,输出

 fun: -0.9999999986126238nit: 2x: array([-5.4602e-08,  1.0000e+00, -8.5667e-09])

意为phr只用了2次迭代,即得到最优解 x 0 = ( 0 1 0 ) \boldsymbol{x}_0=\begin{pmatrix}0\\1\\0\end{pmatrix} x0= 010 ,最优值 f ( x 0 ) = − 1 f(\boldsymbol{x}_0)=-1 f(x0)=1。而对同一问题,sumt(见博文《最优化方法Python计算:求解约束优化问题的罚函数算法》)却用了8次迭代获得同样的结果。
例3:用phr函数求解二次规划
{ minimize x 1 2 + x 1 x 2 + 2 x 2 2 + x 3 2 − 6 x 1 − 2 x 2 − 12 x 3 s.t.   x 1 + x 2 + x 3 − 2 = 0 x 1 − 2 x 2 + 3 ≥ 0 x 1 , x 2 , x 3 ≥ 0 , \begin{cases} \text{minimize}\quad x_1^2+x_1x_2+2x_2^2+x_3^2-6x_1-2x_2-12x_3\\ \text{s.t.\ \ }\quad\quad\quad x_1+x_2+x_3-2=0\\ \quad\quad\quad\quad\quad x_1-2x_2+3\geq0\\ \quad\quad\quad\quad\quad x_1,x_2,x_3\geq0 \end{cases}, minimizex12+x1x2+2x22+x326x12x212x3s.t.  x1+x2+x32=0x12x2+30x1,x2,x30,
给定初始可行解 x 1 = ( 1 1 0 ) \boldsymbol{x}_1=\begin{pmatrix}1\\1\\0\end{pmatrix} x1= 110
:本二次规划的数据矩阵为
H = ( 2 1 0 1 4 0 0 0 2 ) , c = ( − 6 − 2 − 12 ) A e q = ( 1 1 1 ) , b e q = ( 2 ) A i q = ( 1 − 2 0 1 0 0 0 1 0 0 0 1 ) , b i q = ( − 3 0 0 0 ) \begin{array}{l} \boldsymbol{H}=\begin{pmatrix}2&1&0\\1&4&0\\0&0&2\end{pmatrix},\boldsymbol{c}=\begin{pmatrix}-6\\-2\\-12\end{pmatrix}\\ \boldsymbol{A}_{eq}=\begin{pmatrix}1&1&1\end{pmatrix},\boldsymbol{b}_{eq}=(2)\\ \boldsymbol{A}_{iq}=\begin{pmatrix}1&-2&0\\1&0&0\\0&1&0\\0&0&1\end{pmatrix},\boldsymbol{b}_{iq}=\begin{pmatrix}-3\\0\\0\\0\end{pmatrix} \end{array} H= 210140002 ,c= 6212 Aeq=(111),beq=(2)Aiq= 110020100001 ,biq= 3000
下列代码完成计算。

import numpy as np									#导入numpy
H = np.array([[2, 1, 0],							#矩阵H[1, 4, 0],[0, 0, 2]])
c = np.array([-6, -2, -12])							#向量c
Ae = np.array([[1, 1, 1]])							#矩阵Aeq
be = np.array([2])									#向量beq
Ai = np.array([[-1, -2, 0],							#矩阵Aiq[1, 0, 0],[0, 1, 0],[0, 0, 1]])
bi = np.array([-3, 0, 0, 0])						#向量biq
x1 = np.array([1, 1, 0])							#初始迭代点
f = lambda x: 0.5 * np.matmul(x, np.matmul(H, x)) + np.matmul(c, x)	#目标函数
h = lambda x: np.matmul(Ae, x) - be					#等式约束函数
g = lambda x: np.matmul(Ai, x) - bi					#不等式约束函数
print(phr(f, x1, h, g))
\end{lstlisting}\par

运行程序,输出

 fun: -20.000001376228557nit: 7x: array([-9.6384e-08, -1.2540e-07,  2.0000e+00])

即phr函数用7次迭代,算得该二次规划的最优解 x 0 = ( 0 0 2 ) \boldsymbol{x}_0=\begin{pmatrix}0\\0\\2\end{pmatrix} x0= 002 ,最优值 f ( x 0 ) = − 20 f(\boldsymbol{x}_0)=-20 f(x0)=20。相较于sumt函数(见博文《最优化方法Python计算:求解约束优化问题的罚函数算法》)用16次迭代获得同样结果,效率提高是明显的。
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/53289.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Jenkins+docker+springboot 一键自动部署项目步骤

Jenkins+docker+springboot 一键自动部署项目步骤 部署环境:centos+gitee+docker+springboot 实现步骤: docker 安装 jenkins配置 jenkins利用 Dockerfile 和 shell 脚本实现项目打包运行安装 docker docker 安装社区版本 CE 安装需要的软件包yum install -y yum-utils de…

模型 涌现思想

系列文章 分享 模型&#xff0c;了解更多&#x1f449; 模型_思维模型目录。整体产生新特性&#xff0c;超越部分之和。 1 涌现思想的应用 1.1 蚁群算法中的涌现思想 蚁群算法&#xff08;Ant Colony Optimization, ACO&#xff09;是一种模拟蚂蚁觅食行为的计算模型&#xf…

QT项目实战之音乐播放器2.0版本

该版本相较于1.0版本最主要的不同在于连接数据库实现类似于歌曲收藏和取消收藏的功能。 详细情况看我的这篇文章http://t.csdnimg.cn/WS5s8。 效果展示 VSMyMusicShow2.0 define.h UseMySQL.h musicInfo.h VSMyMusicPlayer.h

PMP–一、二、三模–分类–14.敏捷–技巧–帮助团队交付价值的执行实践迭代和增量如何帮助交付工作产品

文章目录 技巧一模14.敏捷--实践--帮助团队交付价值的执行实践--持续集成--在不同层面测试、验收测试驱动开发 (ATDD) 、测试驱动开发和行为驱动开发、刺探 。90、 [单选] 敏捷项目的第一次迭代即将开始。发起人召集团队、Scrum主管、产品负责人和其他项目干系人参加启动会议。…

木舟0基础学习Java的第二十六天(JavaWeb)

设置响应头 resp.setHeader("key","nihao");//推荐使用英文 中文会乱码 案例&#xff1a;模拟登录 jdbc.properties driverClasscom.mysql.jdbc.Driver urljdbc:mysql://localhost:3306/test?verifyServerCertificatefalse&useSSLfalse nameroot p…

第十三届山东省ICPC

vp链接&#xff1a;https://codeforces.com/gym/104417 A. Orders 根据题意模拟&#xff0c;分别按照 a&#xff0c;b 排序&#xff0c;排序后再判断该订单是否能完成。 #include <bits/stdc.h> using namespace std;#define int long longconst int N 105; int n, k…

pikachu文件包含漏洞靶场

本地文件包含 1、先随意进行提交 可以得出是GET传参 可以在filename参数进行文件包含 2、准备一个2.jpg文件 内容为<?php phpinfo();?> 3、上传2.jpg文件 4、访问文件保存的路径uploads/2.jpg 5、将我们上传的文件包含进来 使用../返回上级目录 来进行包含木马文件 …

备战秋招60天算法挑战,Day33

题目链接&#xff1a; https://leetcode.cn/problems/longest-increasing-subsequence/ 视频题解&#xff1a; https://www.bilibili.com/video/BV1RRvheFEog/ LeetCode 300. 最长递增子序列 题目描述 给你一个整数数组nums &#xff0c;找到其中最长严格递增子序列的长度。 …

自幂数判断c++

题目描述 样例输入 3 152 111 153样例输出 F F T 代码如下&#xff1a; #include<bits/stdc.h> using namespace std; long long m,a; int main(){cin>>m;for(int i1;i<m;i){cin>>a;long long ta,n[10],cc0,s0;while(t!0){//求位数与拆位n[cc]t%10;tt/…

多线程篇(并发相关类- 原子操作类)(持续更新迭代)

目录 前言 一、原子变量操作类&#xff08;AtomicLong为例&#xff09; 1. 前言 2. 实例 二、JDK 8新增的原子操作类LongAdder 三、LongAccumulator类原理探究 前言 JUC包提供了一系列的原子性操作类&#xff0c;这些类都是使用非阻塞算法CAS实现的&#xff0c;相比使用…

PowerShell脚本编写:自动化Windows开发工作流程实例介绍

PowerShell脚本编写&#xff1a;自动化Windows开发工作流程 大家好&#xff0c;我是吃点李子&#xff0c;今天我想和大家分享一下如何利用PowerShell脚本编写来自动化我们的日常工作流程。作为一名开发者&#xff0c;我们经常会遇到一些重复性的工作&#xff0c;比如部署环境、…

一起学习LeetCode热题100道(72/100)

72.每日温度(学习) 给定一个整数数组 temperatures &#xff0c;表示每天的温度&#xff0c;返回一个数组 answer &#xff0c;其中 answer[i] 是指对于第 i 天&#xff0c;下一个更高温度出现在几天后。如果气温在这之后都不会升高&#xff0c;请在该位置用 0 来代替。 示例…

【每日一题】LeetCode 2379.得到K个黑块的最少涂色次数(字符串、滑动窗口)

【每日一题】LeetCode 2379.得到K个黑块的最少涂色次数&#xff08;字符串、滑动窗口&#xff09; 题目描述 给定一个字符串 blocks&#xff0c;其中每个字符代表一个颜色块&#xff0c;可以是 ‘W’&#xff08;白色&#xff09;或 ‘B’&#xff08;黑色&#xff09;。你需…

dubbo 服务消费原理分析之应用级服务发现

文章目录 前言一、MigrationRuleListener1、迁移状态模型2、Provider 端升级3、Consumer 端升级4、服务消费选址5、MigrationRuleListener.onRefer6、MigrationRuleHandler.doMigrate6、MigrationRuleHandler.refreshInvoker7、MigrationClusterInvoker.migrateToApplicationFi…

硬件工程师求职过程的启示

在技术求职过程中&#xff0c;每一次面试都是一次独特的学习经历&#xff0c;尤其是在硬件工程师的岗位中&#xff0c;面试者往往需要具备广泛的技能和灵活的适应能力。 1. 面试中的全面技能要求 在现代技术行业中&#xff0c;特别是硬件开发岗位&#xff0c;企业对工程师的要…

初识命名空间

1.创建两个命名空间 ip netns add host1 ip netns add host2 2. 查看命名空间 ip netns ls 3 、 创建veth ip -netns host1 link add veth0 type veth peer name host1-peer 4、 查看命名空间接口 ip -netns host1 address 5、 把host1-peer移动到host2命名空间 ip -ne…

ctfshow-nodejs

什么是nodejs Node.js 是一个基于 Chrome V8 引擎的 Javascript 运行环境。可以说nodejs是一个运行环境&#xff0c;或者说是一个 JS 语言解释器 Nodejs 是基于 Chrome 的 V8 引擎开发的一个 C 程序&#xff0c;目的是提供一个 JS 的运行环境。最早 Nodejs 主要是安装在服务器…

SAP B1 基础实操 - 用户定义字段 (UDF)

目录 一、功能介绍 1. 使用场景 2. 操作逻辑 3. 常用定义部分 3.1 主数据 3.2 营销单据 4. 字段设置表单 4.1 字段基础信息 4.2 不同类详细设置 4.3 默认值/必填 二、案例 1 要求 2 操作步骤 一、功能介绍 1. 使用场景 在实施过程中&#xff0c;经常会碰见用户需…

Jmeter使用时小技巧添加“泊松随机定时器“模拟用户思考时间

1、模拟用户思考时间&#xff0c;添加"泊松随机定时器"

HTML 中 a 超链接标签全解析:属性、锚点与伪类

目录 a超链接标签介绍 锚点介绍 页面内锚点 页面外锚点 伪类 a超链接标签介绍 在 HTML 中&#xff0c;<a>标签是一个文本级的标签&#xff0c;同时它还包含伪类&#xff0c;用于根据用户与链接的不同交互状态呈现不同的样式。与之同为文本级标签的<p>&#xff…