最优化方法Python计算:一般凸二次规划的有效集算法

2024-09-05 11:20

本文主要是介绍最优化方法Python计算:一般凸二次规划的有效集算法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

先考虑仅含不等式约束的二次规划
{ minimize 1 2 x ⊤ H x + c ⊤ x s.t.   A x ≥ b . ( 1 ) \begin{cases} \text{minimize}\quad \frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{Ax}\geq\boldsymbol{b} \end{cases}.\quad\quad{(1)} {minimize21xHx+cxs.t.  Axb.(1)
其中, H ∈ R n × n \boldsymbol{H}\in\text{R}^{n\times n} HRn×n对称正定, A = ( a 1 a 2 ⋮ a m ) ∈ R m × n \boldsymbol{A}=\begin{pmatrix}\boldsymbol{a}_1\\\boldsymbol{a}_2\\\vdots\\\boldsymbol{a}_m\end{pmatrix}\in\text{R}^{m\times n} A= a1a2am Rm×n且rank A = m \boldsymbol{A}=m A=m,可行域 Ω = { x ∣ A x = b } \Omega=\{\boldsymbol{x}|\boldsymbol{Ax}=\boldsymbol{b}\} Ω={xAx=b}
解决问题(1)的方法是每次迭代利用有效集将其转换为求解仅含等式约束的二次型问题。具体而言,对当前可行迭代点 x k \boldsymbol{x}_k xk,设在 x k \boldsymbol{x}_k xk处有效集为 I k I_k Ik,即
I k = { i ∣ 1 ≤ i ≤ m , a i x k = b i } . I_k=\{i|1\leq i\leq m,\boldsymbol{a}_i\boldsymbol{x}_k=b_i\}. Ik={i∣1im,aixk=bi}.
f ( x ) = 1 2 x ⊤ H x + c ⊤ x f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x} f(x)=21xHx+cx A k = ( a i ) i ∈ I k \boldsymbol{A}_k=(\boldsymbol{a}_i)_{i\in I_k} Ak=(ai)iIk b k = ( b i ) i ∈ I k \boldsymbol{b}_k=(b_i)_{i\in I_k} bk=(bi)iIk。考虑等式约束二次规划
{ minimize f ( x ) s.t.   A k x = b k . ( 2 ) \begin{cases} \text{minimize}\quad f(\boldsymbol{x})\\ \text{s.t.\ \ }\quad\boldsymbol{A}_k\boldsymbol{x}=\boldsymbol{b}_k \end{cases}.\quad\quad{(2)} {minimizef(x)s.t.  Akx=bk.(2)
称为二次规划(1)的子问题。若设 x = x k + d \boldsymbol{x}=\boldsymbol{x}_k+\boldsymbol{d} x=xk+d,即 d = x − x k \boldsymbol{d}=\boldsymbol{x}-\boldsymbol{x}_k d=xxk。于是
f ( x ) = f ( x k + d ) = 1 2 ( x k + d ) ⊤ H ( x k + d ) + c ⊤ ( x k + d ) = 1 2 d ⊤ H d + x k ⊤ H d + c ⊤ d + 1 2 x k ⊤ H x k + c ⊤ x k = 1 2 d ⊤ H d + d ⊤ ∇ f ( x k ) + f ( x k ) . \begin{align*} f(\boldsymbol{x})&=f(\boldsymbol{x}_k+\boldsymbol{d})=\frac{1}{2}(\boldsymbol{x}_k+\boldsymbol{d})^\top\boldsymbol{H}(\boldsymbol{x}_k+\boldsymbol{d})+\boldsymbol{c}^\top(\boldsymbol{x}_k+\boldsymbol{d})\\ &=\frac{1}{2}\boldsymbol{d}^\top\boldsymbol{Hd}+\boldsymbol{x}_k^\top\boldsymbol{Hd}+\boldsymbol{c}^\top\boldsymbol{d}+\frac{1}{2}\boldsymbol{x}_k^\top\boldsymbol{Hx}_k+\boldsymbol{c}^\top\boldsymbol{x}_k\\ &=\frac{1}{2}\boldsymbol{d}^\top\boldsymbol{Hd}+\boldsymbol{d}^\top\nabla f(\boldsymbol{x}_k)+f(\boldsymbol{x}_k). \end{align*} f(x)=f(xk+d)=21(xk+d)H(xk+d)+c(xk+d)=21dHd+xkHd+cd+21xkHxk+cxk=21dHd+df(xk)+f(xk).
所以,求子问题(2)等价与求解
{ minimize 1 2 d ⊤ H d + ∇ f ( x k ) ⊤ d s.t.   A k d = o . ( 3 ) \begin{cases} \text{minimize}\quad\frac{1}{2}\boldsymbol{d}^\top\boldsymbol{Hd}+\nabla f(\boldsymbol{x}_k)^\top\boldsymbol{d}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{A}_k\boldsymbol{d}=\boldsymbol{o} \end{cases}.\quad\quad{(3)} {minimize21dHd+f(xk)ds.t.  Akd=o.(3)
H \boldsymbol{H} H的正定性及 A k \boldsymbol{A}_k Ak行满秩,可用等式约束二次规划的拉格朗日算法(见博文《最优化方法Python计算:二次规划的拉格朗日算法》)求解问题(3),设最优解
d ∗ = arg ⁡ min ⁡ A k d = o ( 1 2 d ⊤ H d + ∇ f ( x k ) ⊤ d ) \boldsymbol{d}^{*}=\arg\min_{\boldsymbol{A}_k\boldsymbol{d}=\boldsymbol{o}}\left(\frac{1}{2}\boldsymbol{d}^\top\boldsymbol{Hd}+\nabla f(\boldsymbol{x}_k)^\top\boldsymbol{d}\right) d=argAkd=omin(21dHd+f(xk)d)
对应的乘子为 λ ∗ \boldsymbol{\lambda}^{*} λ
(I)若 d ∗ = o \boldsymbol{d}^{*}=\boldsymbol{o} d=o,则 x k \boldsymbol{x}_k xk是子问题(2)的最优解。此时需判断 x k \boldsymbol{x}_k xk是否为原问题(1)的最优解,这可通过检验对应的拉格朗日因子 λ ∗ \boldsymbol{\lambda}^{*} λ是否满足
λ ∗ ≥ o \boldsymbol{\lambda}^{*}\geq\boldsymbol{o} λo
来完成。若上式成立,则 ( x k λ ∗ ) \begin{pmatrix}\boldsymbol{x}_k\\\boldsymbol{\lambda}^{*}\end{pmatrix} (xkλ)为原问题(1)的KKT点。由于 H \boldsymbol{H} H是正定的,故 x k \boldsymbol{x}_k xk即为原问题(1)的最优解。否则,计算,
j = arg ⁡ min ⁡ i ∈ I k { λ i ∗ } , j=\arg\min_{i\in I_k}\{\lambda^{*}_i\}, j=argiIkmin{λi},
剔除有效约束 a j x = a j ( x k + d ) = b j \boldsymbol{a}_j\boldsymbol{x}=\boldsymbol{a}_j(\boldsymbol{x}_k+\boldsymbol{d})=b_j ajx=aj(xk+d)=bj,即使得
{ a j d > 0 a i d = 0 i ∈ I k , i ≠ j . ( ∗ ) \begin{cases} \boldsymbol{a}_j\boldsymbol{d}>0\\ \boldsymbol{a}_i\boldsymbol{d}=0&i\in I_k,i\not=j \end{cases}.\quad\quad{(*)} {ajd>0aid=0iIk,i=j.()
由于 x k \boldsymbol{x}_k xk为子问题(2)的最优解,故
o = ∇ x L ( x k , λ ∗ ) = ∇ f ( x k ) − A k ⊤ λ ∗ \boldsymbol{o}=\nabla_{\boldsymbol{x}}L(\boldsymbol{x}_k,\boldsymbol{\lambda}^{*})=\nabla f(\boldsymbol{x}_k)-\boldsymbol{A}_k^\top\boldsymbol{\lambda}^{*} o=xL(xk,λ)=f(xk)Akλ
其中, L ( x , λ ) = f ( x ) − ( A k x − b k ) ⊤ λ L(\boldsymbol{x},\boldsymbol{\lambda})=f(\boldsymbol{x})-(\boldsymbol{A}_k\boldsymbol{x}-\boldsymbol{b}_k)^\top\boldsymbol{\lambda} L(x,λ)=f(x)(Akxbk)λ为子问题(2)的拉格朗日函数。即
∇ f ( x k ) = A k ⊤ λ ∗ \nabla f(\boldsymbol{x}_k)=\boldsymbol{A}_k^\top\boldsymbol{\lambda}^{*} f(xk)=Akλ
对满足式 ( ∗ ) (*) ()的向量 d \boldsymbol{d} d
∇ f ( x k ) ⊤ d = λ ∗ ⊤ A k d = λ j a j d < 0. \nabla f(\boldsymbol{x}_k)^\top\boldsymbol{d}=\boldsymbol{\lambda}^{*\top}\boldsymbol{A}_k\boldsymbol{d}=\lambda_j\boldsymbol{a}_j\boldsymbol{d}<0. f(xk)d=λAkd=λjajd<0.
d \boldsymbol{d} d是问题(3)的可行下降方向。于是修改有效集下标集
I k = I k − { j } , I_k=I_k-\{j\}, Ik=Ik{j},
重新求解子问题(2),即可得到原问题(1)下降可行方向 d ∗ \boldsymbol{d}^{*} d
d ∗ ≠ o \boldsymbol{d}^{*}\not=\boldsymbol{o} d=o,取搜索方向 d k = d ∗ \boldsymbol{d}_k=\boldsymbol{d}^{*} dk=d
(II)此时,若 x k + d k \boldsymbol{x}_k+\boldsymbol{d}_k xk+dk是原问题(7.4)的可行解,即 x k + d k ∈ Ω \boldsymbol{x}_k+\boldsymbol{d}_k\in\Omega xk+dkΩ,则取
x k + 1 = x k + d k \boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\boldsymbol{d}_k xk+1=xk+dk
作为下一个迭代点。
(III)若 x k + d k ∉ Ω \boldsymbol{x}_k+\boldsymbol{d}_k\not\in\Omega xk+dkΩ,寻求搜索步长 α k \alpha_k αk,使得 x k + α k d k ∈ Ω \boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k\in\Omega xk+αkdkΩ,作为下一个迭代点。为此,需
a i ( x k + α k d k ) ≥ b i , i ∈ I k ‾ . ( ∗ ∗ ) \boldsymbol{a}_i(\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k)\geq b_i,\quad i\in \overline{I_k}.\quad\quad{(**)} ai(xk+αkdk)bi,iIk.()
其中 I k ‾ = { 1 , 2 , ⋯ , m } − I k \overline{I_k}=\{1,2,\cdots,m\}-I_k Ik={1,2,,m}Ik。由于 x k ∈ Ω \boldsymbol{x}_k\in\Omega xkΩ,故 A x k ≥ b k \boldsymbol{Ax}_k\geq\boldsymbol{b}_k Axkbk。因此,若 a i d k ≥ 0 \boldsymbol{a}_i\boldsymbol{d}_k\geq0 aidk0,则对任意 α k > 0 \alpha_k>0 αk>0,式 ( ∗ ∗ ) (**) ()均成立。往设存在 i ∈ I k ‾ i\in \overline{I_k} iIk,使得 a i d k < 0 \boldsymbol{a}_i\boldsymbol{d}_k<0 aidk<0。这时,取
0 < α k ≤ b i − a i x k a i d k 0<\alpha_k\leq\frac{b_i-\boldsymbol{a}_i\boldsymbol{x}_k}{\boldsymbol{a}_i\boldsymbol{d}_k} 0<αkaidkbiaixk
则必有 a i ( x k + α k d k ) ≥ b i \boldsymbol{a}_i(\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k)\geq b_i ai(xk+αkdk)bi。于是,令
α ^ k = min ⁡ ( b i − a i x k a i d k ∣ i ∈ I k ‾ , a i d k < 0 ) \hat{\alpha}_k=\min\left(\frac{b_i-\boldsymbol{a}_i\boldsymbol{x}_k}{\boldsymbol{a}_i\boldsymbol{d}_k}\Big|i\in\overline{I_k},\boldsymbol{a}_i\boldsymbol{d}_k<0\right) α^k=min(aidkbiaixk iIk,aidk<0)
则可保证式
a i ( x k + α ^ k d k ) ≥ b i , i ∈ I k ‾ \boldsymbol{a}_i(\boldsymbol{x}_k+\hat{\alpha}_k\boldsymbol{d}_k)\geq b_i,\quad i\in \overline{I_k} ai(xk+α^kdk)bi,iIk
成立。综合(II)和(III),令
α k = min ⁡ { 1 , α ^ k } \alpha_k=\min\{1,\hat{\alpha}_k\} αk=min{1,α^k}
则当 d k = d ∗ ≠ o \boldsymbol{d}_k=\boldsymbol{d}^{*}\not=\boldsymbol{o} dk=d=o
x k + 1 = x k + α k d k ∈ Ω \boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k\in\Omega xk+1=xk+αkdkΩ
可作为下一个迭代点。此时,需更新 x k + 1 \boldsymbol{x}_{k+1} xk+1的有效集
I k + 1 = { i ∣ 1 ≤ i ≤ m , a i x k + 1 = b i } . I_{k+1}=\{i|1\leq i\leq m,\boldsymbol{a}_i\boldsymbol{x}_{k+1}=b_i\}. Ik+1={i∣1im,aixk+1=bi}.
对既有等式约束亦有不等式约束的二次规划
{ minimize 1 2 x ⊤ H x + c ⊤ x s.t.   A e q x − b e q = o A i q x − b i q ≥ o ( 4 ) \begin{cases} \text{minimize}\quad \frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{A}_{eq}\boldsymbol{x}-\boldsymbol{b}_{eq}=\boldsymbol{o}\\ \quad\quad\quad\quad\quad\boldsymbol{A}_{iq}\boldsymbol{x}-\boldsymbol{b}_{iq}\geq\boldsymbol{o} \end{cases}\quad\quad{(4)} minimize21xHx+cxs.t.  Aeqxbeq=oAiqxbiqo(4)
其中, H \boldsymbol{H} H对称正定, A e q ∈ R l × n \boldsymbol{A}_{eq}\in\text{R}^{l\times n} AeqRl×n A i q ∈ R m × n \boldsymbol{A}_{iq}\in\text{R}^{m\times n} AiqRm×n,且rank ( A e q A i q ) = l + m \begin{pmatrix}\boldsymbol{A}_{eq}\\\boldsymbol{A}_{iq}\end{pmatrix}=l+m (AeqAiq)=l+m。令 A = ( A e q A i q ) = ( a 1 ⋮ a l a l + 1 ⋮ a l + m ) \boldsymbol{A}=\begin{pmatrix}\boldsymbol{A}_{eq}\\\boldsymbol{A}_{iq}\end{pmatrix}=\begin{pmatrix}\boldsymbol{a}_1\\\vdots\\\boldsymbol{a}_l\\\boldsymbol{a}_{l+1}\\\vdots\\\boldsymbol{a}_{l+m}\end{pmatrix} A=(AeqAiq)= a1alal+1al+m b = ( b e q b i q ) = ( b 1 ⋮ b l b l + 1 ⋮ b l + m ) \boldsymbol{b}=\begin{pmatrix}\boldsymbol{b}_{eq}\\\boldsymbol{b}_{iq}\end{pmatrix}=\begin{pmatrix}b_1\\\vdots\\b_l\\b_{l+1}\\\vdots\\b_{l+m}\end{pmatrix} b=(beqbiq)= b1blbl+1bl+m 。记 E = { 1 , ⋯ , l } E=\{1,\cdots,l\} E={1,,l} I = { l + 1 , ⋯ , l + m } I=\{l+1,\cdots,l+m\} I={l+1,,l+m},对任一 x ∈ Ω = { x ∣ A e q x = b e q , A i q x ≥ b i q } \boldsymbol{x}\in\Omega=\{\boldsymbol{x}|\boldsymbol{A}_{eq}\boldsymbol{x}=\boldsymbol{b}_{eq},\boldsymbol{A}_{iq}\boldsymbol{x}\geq\boldsymbol{b}_{iq}\} xΩ={xAeqx=beq,Aiqxbiq} I ( x ) = { i ∣ i ∈ I , a i x = 0 } I(\boldsymbol{x})=\{i|i\in I,\boldsymbol{a}_i\boldsymbol{x}=0\} I(x)={iiI,aix=0} x \boldsymbol{x} x的有效集。 S ( x ) = E ∪ I ( x ) S(\boldsymbol{x})=E\cup I(\boldsymbol{x}) S(x)=EI(x)为子问题(2)的约束下标集。即可用有效约束集方法求解问题(4)。下列代码实现求解一般二次规划(4)的有效集算法。

import numpy as np										#导入numpy
from scipy.optimize import OptimizeResult				#导入OptimizeResult
def qpact(H, c, x1, Aeq = None, beq = None, Aiq = None, biq = None):l = 0												#等式约束个数if isinstance(Aeq, np.ndarray):A = Aeq.copy()b = beq.copy()l, n = Aeq.shapem = 0												#不等式约数个数if isinstance(Aiq, np.ndarray):if isinstance(Aeq, np.ndarray):A = np.vstack((A, Aiq))b = np.concatenate((b, biq))m = biq.sizeelse:A = Aiq.copy()b = biq.copy()m, n = Aiq.shapegk = lambda x: np.matmul(H, x) + c					#目标函数梯度函数E = np.arange(l)									#等式约束下标集I = np.arange(l, l + m)								#不等式约束下标集k = 0												#迭代数xk = x1												#当前迭代点b1 = np.zeros(l + m)Ik = np.where(np.abs(np.matmul(A[l:,:],xk) - b[l:])<1e-10)[0] + l	#有效集Sk = np.concatenate((E, Ik))						#子问题约束opt = Falsewhile not opt:dk, lamdk = qlag(H, A[Sk], b1[Sk], gk(xk))		#计算搜索方向if np.linalg.norm(dk) > 1e-6:					#搜索方向非零Ikb = np.setdiff1d(I, Ik)					#非有效集neg = np.where(np.matmul(A[Ikb], dk) < 0)[0]arr = (b[Ikb[neg]] - np.matmul(A[Ikb[neg]], xk)) / np.matmul(A[Ikb[neg]], dk)alpha1 = 1if arr.size > 0:alpha1 = np.min(arr)alpha = min(1, alpha1)						#搜索步长xk = xk + alpha * dkIk = np.where(np.abs(np.matmul(A[l:,:], xk) - b[l:]) < 1e-10)[0] + l#重置有效集Sk = np.concatenate((E, Ik))				#重置子问题约束k += 1else:minlamd = 0.0if lamdk.size > l:minlamd = np.min(lamdk[l:])if minlamd < 0:								#乘子存在负值元素j = np.argmin(lamdk[l:])Ik = np.delete(Ik, j)					#修改有效集Sk = np.concatenate((E, Ik))			#更新子问题约束集else:										#乘子非负opt = Truebestx = xkbesty = 0.5 * np.matmul(xk, np.matmul(H, xk)) + np.matmul(c, xk)return OptimizeResult(fun = besty, x = bestx, nit = k)

程序的第3~54行定义的qpact函数实现求解二次规划问题(4)的有效集算法。该函数的参数H,c,Aeq,beq,Aiq,biq分别表示问题(1)
{ minimize 1 2 x ⊤ H x + c ⊤ x s.t.   A e q x − b e q = o A i q x − b i q ≥ o \begin{cases} \text{minimize}\quad \frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{A}_{eq}\boldsymbol{x}-\boldsymbol{b}_{eq}=\boldsymbol{o}\\ \quad\quad\quad\quad\quad\boldsymbol{A}_{iq}\boldsymbol{x}-\boldsymbol{b}_{iq}\geq\boldsymbol{o} \end{cases} minimize21xHx+cxs.t.  Aeqxbeq=oAiqxbiqo
中的 H \boldsymbol{H} H c \boldsymbol{c} c A e q \boldsymbol{A}_{eq} Aeq b e q \boldsymbol{b}_{eq} beq A i q \boldsymbol{A}_{iq} Aiq b i q \boldsymbol{b}_{iq} biq,x1表示问题的初始可行解 x 1 \boldsymbol{x}_1 x1。其中,Aeq,beq,Aiq和biq是命名参数,缺省值为无类型常量None。
函数体内,第4~{}18行将表示 A e q \boldsymbol{A}_{eq} Aeq A i q \boldsymbol{A}_{iq} Aiq的Aeq和Aiq叠加成 A = ( A e q A i q ) \boldsymbol{A}=\begin{pmatrix}\boldsymbol{A}_{eq}\\\boldsymbol{A}_{iq}\end{pmatrix} A=(AeqAiq),赋予A。将表示 b e q \boldsymbol{b}_{eq} beq b i q \boldsymbol{b}_{iq} biq的beq和biq叠加成 b = ( b e q b i q ) \boldsymbol{b}=\begin{pmatrix}\boldsymbol{b}_{eq}\\\boldsymbol{b}_{iq}\end{pmatrix} b=(beqbiq),赋予b。记录下等式约束个数l和不等式约束个数m。
第19行定义问题(4)的目标函数的梯度 ∇ f ( x ) = H x + c \nabla f(\boldsymbol{x})=\boldsymbol{Hx}+\boldsymbol{c} f(x)=Hx+c赋予gk。第20、21行设置等式约束下标集E和不等式约束下标集I。第22~24初始化迭代次数k,迭代点xk,零向量b1。第25初始化xk处的有效集Ik,第26行初始化子问题的等式约束下标集Sk。
第28~51行的while循环,执行算法的迭代运算。其中,第29行调用博文《最优化方法Python计算:二次规划的拉格朗日算法》定义的拉格朗日算法函数qlag,求解
{ minimize 1 2 d ⊤ H d + ∇ f ( x k ) ⊤ d s.t.   a i d = 0 , i ∈ I k , \begin{cases} \text{minimize}\quad\frac{1}{2}\boldsymbol{d}^\top\boldsymbol{Hd}+\nabla f(\boldsymbol{x}_k)^\top\boldsymbol{d}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{a}_i\boldsymbol{d}=0,\quad i\in I_k \end{cases}, {minimize21dHd+f(xk)ds.t.  aid=0,iIk,
以求搜索方向 d k \boldsymbol{d}_k dk及其对应的拉格朗日乘子 λ k \boldsymbol{\lambda}_k λk,赋予dk和lamdk。第30~41行处理当搜索方向 d k ≠ o \boldsymbol{d}_k\not=\boldsymbol{o} dk=o时,按
α k = min ⁡ { 1 , α ^ k } \alpha_k=\min\{1,\hat{\alpha}_k\} αk=min{1,α^k}
其中, α ^ k = min ⁡ ( b i − a i x k a i d k ∣ i ∈ I k ‾ , a i d k < 0 ) \hat{\alpha}_k=\min\left(\frac{b_i-\boldsymbol{a}_i\boldsymbol{x}_k}{\boldsymbol{a}_i\boldsymbol{d}_k}\Big|i\in\overline{I_k},\boldsymbol{a}_i\boldsymbol{d}_k<0\right) α^k=min(aidkbiaixk iIk,aidk<0),计算搜索步长 α k \alpha_k αk,赋予alpha。以xk+alpha*dk更新迭代点xk,并更新有效集Ik和子问题等式约束下标集Sk。第42~51行处理所得搜索方向 d k = o \boldsymbol{d}_k=\boldsymbol{o} dk=o的情形。根据拉格朗日乘子 λ k \boldsymbol{\lambda}_k λk是否存在负值元素,或更新有效集并进入下一轮迭代(第44~49行)或完成迭代置标志opt为True(第51行)。第52~54行返回最优解xk及最优目标近似值。
例1: 用qpact函数求二次规划
{ minimize x 1 2 + 2 x 2 2 + x 3 2 − 2 x 1 x 2 + x 3 s.t.   x 1 + x 2 + x 3 = 4 2 x 1 − x 2 + x 3 = 2 , \begin{cases} \text{minimize}\quad x_1^2+2x_2^2+x_3^2-2x_1x_2+x_3\\ \text{s.t.\ \ }\quad\quad\quad x_1+x_2+x_3=4\\ \quad\quad\quad\quad\quad 2x_1-x_2+x_3=2 \end{cases}, minimizex12+2x22+x322x1x2+x3s.t.  x1+x2+x3=42x1x2+x3=2,
给定初始可行解 x 1 = ( 2 2 0 ) \boldsymbol{x}_1=\begin{pmatrix}2\\2\\0\end{pmatrix} x1= 220
:本题中,
H = ( 2 − 2 0 − 2 4 0 0 0 2 ) , A e q = ( 1 1 1 2 − 1 1 ) , b e q = ( 4 2 ) , c = ( 0 0 1 ) \boldsymbol{H}=\begin{pmatrix}2&-2&0\\-2&4&0\\0&0&2\end{pmatrix},\boldsymbol{A}_{eq}=\begin{pmatrix}1&1&1\\2&-1&1\end{pmatrix},\boldsymbol{b}_{eq}=\begin{pmatrix}4\\2\end{pmatrix},\boldsymbol{c}=\begin{pmatrix}0\\0\\1\end{pmatrix} H= 220240002 ,Aeq=(121111),beq=(42),c= 001
下列代码完成计算。

import numpy as np						#导入numpy
from fractions import Fraction as F		#设置输出格式
np.set_printoptions(formatter={'all':lambda x:str(F(x).limit_denominator())})
H = np.array([[2, -2, 0],				#矩阵H[-2, 4, 0],[0, 0, 2]])
c = np.array([0, 0, 1])					#向量c
Ae = np.array([[1, 1, 1],				#矩阵Aeq[2, -1, 1]])
be = np.array([4, 2])					#向量beq
x = np.array([2, 2, 0])					#初始可行解
print(qpact(H, c, x, Ae, be))

借助代码内部注释信息,不难理解本程序。需要注意的是本问题只有等式约束,没有不等式约束。所以,第14行调用qpact函数时向形参Aeq和beq(位于形参Aiq和biq之前)直接传递实际参数Ae、be即可。参数Aiq和biq使用缺省值。运行程序,输出

 fun: 3.977272727272721nit: 1x: array([21/11, 43/22, 3/22])

仅迭代1次即算得最优解 x 0 = ( 21 11 43 22 3 22 ) \boldsymbol{x}_0=\begin{pmatrix}\frac{21}{11}\\\frac{43}{22}\\\frac{3}{22}\end{pmatrix} x0= 11212243223
例2:用qpact函数求解二次规划
{ minimize x 1 2 − x 1 x 2 + 2 x 2 2 − x 1 − 10 x 2 s.t.   − 3 x 1 − 2 x 2 ≥ − 6 x 1 , x 2 ≥ 0 , \begin{cases} \text{minimize}\quad x_1^2-x_1x_2+2x_2^2-x_1-10x_2\\ \text{s.t.\ \ }\quad\quad -3x_1-2x_2\geq-6\\ \quad\quad\quad\quad\quad x_1,x_2\geq0 \end{cases}, minimizex12x1x2+2x22x110x2s.t.  3x12x26x1,x20,
给定初始可行解 x 1 = o \boldsymbol{x}_1=\boldsymbol{o} x1=o
:本题中
H = ( 2 − 1 − 1 4 ) , c = ( − 1 − 10 ) , x 1 = ( 0 0 ) , A i q = ( − 3 − 2 1 0 0 1 ) , b i q = ( − 6 0 0 ) . \boldsymbol{H}=\begin{pmatrix}2&-1\\-1&4\end{pmatrix},\boldsymbol{c}=\begin{pmatrix}-1\\-10\end{pmatrix},\boldsymbol{x}_1=\begin{pmatrix}0\\0\end{pmatrix},\boldsymbol{A}_{iq}=\begin{pmatrix}-3&-2\\1&0\\0&1\end{pmatrix},\boldsymbol{b}_{iq}=\begin{pmatrix}-6\\0\\0\end{pmatrix}. H=(2114),c=(110),x1=(00),Aiq= 310201 ,biq= 600 .
下列代码利用这些数据完成本例计算。

import numpy as np					#导入numpy
from utility import qpact			#导入qpact
from fractions import Fraction as F	#设置输出格式
np.set_printoptions(formatter={'all':lambda x:str(F(x).limit_denominator())})
H = np.array([[2, -1],				#矩阵H[-1, 4]])
c = np.array([-1, -10])				#向量c
Ai = np.array([[-3, -2],			#矩阵Ai[1, 0],[0, 1]])
bi = np.array([-6, 0, 0])			#向量bi
x = np.array([0, 0])				#初始可行解
print(qpact(H, c, x, Aiq=Ai, biq=bi))

注意本问题只有不等式约束,没有等式约束。所以,第13行调用qpact函数传递实际参数Ai、bi时,须用赋值形式将其传递给形式参数Aiq和biq(位于形参Aeq和beq之后)。参数Aeq和beq使用缺省值。运行程序,输出

 fun: -13.75nit: 3x: array([1/2, 9/4])

这意味着,经过3此迭代,有效集算法算得本问题最优解 x 0 = ( 1 2 9 4 ) \boldsymbol{x}_0=\begin{pmatrix}\frac{1}{2}\\\frac{9}{4}\end{pmatrix} x0=(2149),最优值 f ( x 0 ) ≈ − 13.75 f(\boldsymbol{x}_0)\approx-13.75 f(x0)13.75
例3:用qpact函数求解二次规划
{ 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
{\heiti{解}}:本问题中的数据矩阵表示为
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
from utility import qpact				#导入qpact
from fractions import Fraction as F		#设置输出格式
np.set_printoptions(formatter={'all':lambda x:str(F(x).limit_denominator())})
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向量
x = np.array([1, 1, 0])					#初始可行解
print(qpact(H, c, x, Ae, be, Ai, bi))

由于本问题既有等式约束,又有不等式约束,第17行调用qpact时,所有实际参数按形式参数的顺序传递即可。运行程序,输出

 fun: -20.0nit: 2x: array([0, 0, 2])

即2次迭代,算得最优解为 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
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!

这篇关于最优化方法Python计算:一般凸二次规划的有效集算法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1138818

相关文章

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后,更是加大了理解的难度。 但,不了解推荐算法,就无法做推荐系

python: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别: 当你在一个模块(文件)中使用 global 声明变量时,这个变量只在该模块的全局命名空

康拓展开(hash算法中会用到)

康拓展开是一个全排列到一个自然数的双射(也就是某个全排列与某个自然数一一对应) 公式: X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。(a[i]在不同应用中的含义不同); 典型应用: 计算当前排列在所有由小到大全排列中的顺序,也就是说求当前排列是第

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动,这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升:减少抖动,提高图像的清晰度和细节表现力,使得监控画面更加真实可信。 - 细节增强:在低光条件下,抖

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

动态规划---打家劫舍

题目: 你是一个专业的小偷,计划偷窃沿街的房屋。每间房内都藏有一定的现金,影响你偷窃的唯一制约因素就是相邻的房屋装有相互连通的防盗系统,如果两间相邻的房屋在同一晚上被小偷闯入,系统会自动报警。 给定一个代表每个房屋存放金额的非负整数数组,计算你 不触动警报装置的情况下 ,一夜之内能够偷窃到的最高金额。 思路: 动态规划五部曲: 1.确定dp数组及含义 dp数组是一维数组,dp[i]代表

【C++】_list常用方法解析及模拟实现

相信自己的力量,只要对自己始终保持信心,尽自己最大努力去完成任何事,就算事情最终结果是失败了,努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一:什么是list? •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二:list常用接口 •🌰1.默认成员函数 🔥构造函数(⭐) 🔥析构函数 •🌰2.list对象

【Python编程】Linux创建虚拟环境并配置与notebook相连接

1.创建 使用 venv 创建虚拟环境。例如,在当前目录下创建一个名为 myenv 的虚拟环境: python3 -m venv myenv 2.激活 激活虚拟环境使其成为当前终端会话的活动环境。运行: source myenv/bin/activate 3.与notebook连接 在虚拟环境中,使用 pip 安装 Jupyter 和 ipykernel: pip instal

浅谈主机加固,六种有效的主机加固方法

在数字化时代,数据的价值不言而喻,但随之而来的安全威胁也日益严峻。从勒索病毒到内部泄露,企业的数据安全面临着前所未有的挑战。为了应对这些挑战,一种全新的主机加固解决方案应运而生。 MCK主机加固解决方案,采用先进的安全容器中间件技术,构建起一套内核级的纵深立体防护体系。这一体系突破了传统安全防护的局限,即使在管理员权限被恶意利用的情况下,也能确保服务器的安全稳定运行。 普适主机加固措施: