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

先考虑仅含不等式约束的二次规划
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }\mathbf{Ax}\ge \mathbf{b}\end{array}.(1)$$
其中，$\mathbf{H}\in {\text{R}}^{n\times n}$对称正定，$\mathbf{A}=\left(\begin{array}{c}{\mathbf{a}}_1\\ {\mathbf{a}}_2\\ ⋮\\ {\mathbf{a}}_m\end{array}\right)\in {\text{R}}^{m\times n}$且rank$\mathbf{A}=m$，可行域 Ω = { x ∣ A x = b } \Omega=\{\boldsymbol{x}|\boldsymbol{Ax}=\boldsymbol{b}\} Ω={x∣Ax=b}。
解决问题(1)的方法是每次迭代利用有效集将其转换为求解仅含等式约束的二次型问题。具体而言，对当前可行迭代点${\mathbf{x}}_k$，设在${\mathbf{x}}_k$处有效集为$I_k$，即
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∣1≤i≤m,ai​xk​=bi​}.
记$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}$，${\mathbf{A}}_k=({\mathbf{a}}_i{)}_{i\in I_k}$，${\mathbf{b}}_k=(b_i{)}_{i\in I_k}$。考虑等式约束二次规划
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }{\mathbf{A}}_k\mathbf{x}={\mathbf{b}}_k\end{array}.(2)$$
称为二次规划(1)的子问题。若设$\mathbf{x}={\mathbf{x}}_k+\mathbf{d}$，即$\mathbf{d}=\mathbf{x}-{\mathbf{x}}_k$。于是
$$\begin{array}{rl}f(\mathbf{x})& =f({\mathbf{x}}_k+\mathbf{d})=\frac{1}{2}({\mathbf{x}}_k+\mathbf{d}{)}^{\top }\mathbf{H}({\mathbf{x}}_k+\mathbf{d})+{\mathbf{c}}^{\top }({\mathbf{x}}_k+\mathbf{d})\\ & =\frac{1}{2}{\mathbf{d}}^{\top }\mathbf{Hd}+{\mathbf{x}}_k^{\top }\mathbf{Hd}+{\mathbf{c}}^{\top }\mathbf{d}+\frac{1}{2}{\mathbf{x}}_k^{\top }{\mathbf{Hx}}_k+{\mathbf{c}}^{\top }{\mathbf{x}}_k\\ & =\frac{1}{2}{\mathbf{d}}^{\top }\mathbf{Hd}+{\mathbf{d}}^{\top }\nabla f({\mathbf{x}}_k)+f({\mathbf{x}}_k).\end{array}$$
所以，求子问题(2)等价与求解
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{d}}^{\top }\mathbf{Hd}+\nabla f({\mathbf{x}}_k{)}^{\top }\mathbf{d}\\ \text{s.t. }{\mathbf{A}}_k\mathbf{d}=\mathbf{o}\end{array}.(3)$$
由$\mathbf{H}$的正定性及${\mathbf{A}}_k$行满秩，可用等式约束二次规划的拉格朗日算法（见博文《最优化方法Python计算：二次规划的拉格朗日算法》）求解问题(3)，设最优解
$${\mathbf{d}}^{\ast }=\arg \underset{{\mathbf{A}}_k\mathbf{d}=\mathbf{o}}{\min }\left(\frac{1}{2}{\mathbf{d}}^{\top }\mathbf{Hd}+\nabla f({\mathbf{x}}_k{)}^{\top }\mathbf{d}\right)$$
对应的乘子为${\mathbf{\lambda }}^{\ast }$。
(I)若${\mathbf{d}}^{\ast }=\mathbf{o}$，则${\mathbf{x}}_k$是子问题(2)的最优解。此时需判断${\mathbf{x}}_k$是否为原问题(1)的最优解，这可通过检验对应的拉格朗日因子${\mathbf{\lambda }}^{\ast }$是否满足
$${\mathbf{\lambda }}^{\ast }\ge \mathbf{o}$$
来完成。若上式成立，则$\left(\begin{array}{c}{\mathbf{x}}_k\\ {\mathbf{\lambda }}^{\ast }\end{array}\right)$为原问题(1)的KKT点。由于$\mathbf{H}$是正定的，故${\mathbf{x}}_k$即为原问题(1)的最优解。否则，计算，
j = arg ⁡ min ⁡ i ∈ I k { λ i ∗ } , j=\arg\min_{i\in I_k}\{\lambda^{*}_i\}, j=argi∈Ik​min​{λi∗​},
剔除有效约束${\mathbf{a}}_j\mathbf{x}={\mathbf{a}}_j({\mathbf{x}}_k+\mathbf{d})=b_j$，即使得
$$\{\begin{array}{ll}{\mathbf{a}}_j\mathbf{d}>0& \\ {\mathbf{a}}_i\mathbf{d}=0& i\in I_k,i\ne j\end{array}.(\ast )$$
由于${\mathbf{x}}_k$为子问题(2)的最优解，故
$$\mathbf{o}={\nabla }_{\mathbf{x}}L({\mathbf{x}}_k,{\mathbf{\lambda }}^{\ast })=\nabla f({\mathbf{x}}_k)-{\mathbf{A}}_k^{\top }{\mathbf{\lambda }}^{\ast }$$
其中，$L(\mathbf{x},\mathbf{\lambda })=f(\mathbf{x})-({\mathbf{A}}_k\mathbf{x}-{\mathbf{b}}_k{)}^{\top }\mathbf{\lambda }$为子问题(2)的拉格朗日函数。即
$$\nabla f({\mathbf{x}}_k)={\mathbf{A}}_k^{\top }{\mathbf{\lambda }}^{\ast }$$
对满足式$(\ast )$的向量$\mathbf{d}$，
$$\nabla f({\mathbf{x}}_k{)}^{\top }\mathbf{d}={\mathbf{\lambda }}^{\ast \top }{\mathbf{A}}_k\mathbf{d}={\lambda }_j{\mathbf{a}}_j\mathbf{d}<0.$$
即$\mathbf{d}$是问题(3)的可行下降方向。于是修改有效集下标集
I k = I k − { j } , I_k=I_k-\{j\}, Ik​=Ik​−{j},
重新求解子问题(2)，即可得到原问题(1)下降可行方向${\mathbf{d}}^{\ast }$。
若${\mathbf{d}}^{\ast }\ne \mathbf{o}$，取搜索方向${\mathbf{d}}_k={\mathbf{d}}^{\ast }$。
(II)此时，若${\mathbf{x}}_k+{\mathbf{d}}_k$是原问题(7.4)的可行解，即${\mathbf{x}}_k+{\mathbf{d}}_k\in \Omega$，则取
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{d}}_k$$
作为下一个迭代点。
(III)若${\mathbf{x}}_k+{\mathbf{d}}_k\notin \Omega$，寻求搜索步长${\alpha }_k$，使得${\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k\in \Omega$，作为下一个迭代点。为此，需
$${\mathbf{a}}_i({\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k)\ge b_i,i\in \overline{I_k}.(\ast \ast )$$
其中 I k ‾ = { 1 , 2 , ⋯   , m } − I k \overline{I_k}=\{1,2,\cdots,m\}-I_k Ik​​={1,2,⋯,m}−Ik​。由于${\mathbf{x}}_k\in \Omega$，故${\mathbf{Ax}}_k\ge {\mathbf{b}}_k$。因此，若${\mathbf{a}}_i{\mathbf{d}}_k\ge 0$，则对任意${\alpha }_k>0$，式$(\ast \ast )$均成立。往设存在$i\in \overline{I_k}$，使得${\mathbf{a}}_i{\mathbf{d}}_k<0$。这时，取
$$0<{\alpha }_k\le \frac{b_i-{\mathbf{a}}_i{\mathbf{x}}_k}{{\mathbf{a}}_i{\mathbf{d}}_k}$$
则必有${\mathbf{a}}_i({\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k)\ge b_i$。于是，令
$${\hat{\alpha }}_k=\min \left(\frac{b_i-{\mathbf{a}}_i{\mathbf{x}}_k}{{\mathbf{a}}_i{\mathbf{d}}_k}|i\in \overline{I_k},{\mathbf{a}}_i{\mathbf{d}}_k<0\right)$$
则可保证式
$${\mathbf{a}}_i({\mathbf{x}}_k+{\hat{\alpha }}_k{\mathbf{d}}_k)\ge b_i,i\in \overline{I_k}$$
成立。综合(II)和(III)，令
α k = min ⁡ { 1 , α ^ k } \alpha_k=\min\{1,\hat{\alpha}_k\} αk​=min{1,α^k​}
则当${\mathbf{d}}_k={\mathbf{d}}^{\ast }\ne \mathbf{o}$时
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k\in \Omega$$
可作为下一个迭代点。此时，需更新${\mathbf{x}}_{k+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∣1≤i≤m,ai​xk+1​=bi​}.
对既有等式约束亦有不等式约束的二次规划
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }{\mathbf{A}}_{eq}\mathbf{x}-{\mathbf{b}}_{eq}=\mathbf{o}\\ {\mathbf{A}}_{iq}\mathbf{x}-{\mathbf{b}}_{iq}\ge \mathbf{o}\end{array}(4)$$
其中，$\mathbf{H}$对称正定，${\mathbf{A}}_{eq}\in {\text{R}}^{l\times n}$，${\mathbf{A}}_{iq}\in {\text{R}}^{m\times n}$，且rank$\left(\begin{array}{c}{\mathbf{A}}_{eq}\\ {\mathbf{A}}_{iq}\end{array}\right)=l+m$。令$\mathbf{A}=\left(\begin{array}{c}{\mathbf{A}}_{eq}\\ {\mathbf{A}}_{iq}\end{array}\right)=\left(\begin{array}{c}{\mathbf{a}}_1\\ ⋮\\ {\mathbf{a}}_l\\ {\mathbf{a}}_{l+1}\\ ⋮\\ {\mathbf{a}}_{l+m}\end{array}\right)$，$\mathbf{b}=\left(\begin{array}{c}{\mathbf{b}}_{eq}\\ {\mathbf{b}}_{iq}\end{array}\right)=\left(\begin{array}{c}{b}_1\\ ⋮\\ b_l\\ b_{l+1}\\ ⋮\\ b_{l+m}\end{array}\right)$。记 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∈Ω={x∣Aeq​x=beq​,Aiq​x≥biq​}， I ( x ) = { i ∣ i ∈ I , a i x = 0 } I(\boldsymbol{x})=\{i|i\in I,\boldsymbol{a}_i\boldsymbol{x}=0\} I(x)={i∣i∈I,ai​x=0}为$\mathbf{x}$的有效集。$S(\mathbf{x})=E\cup I(\mathbf{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.shape
    m = 0												#不等式约数个数
    if isinstance(Aiq, np.ndarray):
        if isinstance(Aeq, np.ndarray):
            A = np.vstack((A, Aiq))
            b = np.concatenate((b, biq))
            m = biq.size
        else:
            A = Aiq.copy()
            b = biq.copy()
            m, n = Aiq.shape
    gk = 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 = False
    while 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 = 1
            if arr.size > 0:
                alpha1 = np.min(arr)
            alpha = min(1, alpha1)						#搜索步长
            xk = xk + alpha * dk
            Ik = np.where(np.abs(np.matmul(A[l:,:], xk) - b[l:]) < 1e-10)[0] + l#重置有效集
            Sk = np.concatenate((E, Ik))				#重置子问题约束
            k += 1
        else:
            minlamd = 0.0
            if 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 = True
    bestx = xk
    besty = 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)
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }{\mathbf{A}}_{eq}\mathbf{x}-{\mathbf{b}}_{eq}=\mathbf{o}\\ {\mathbf{A}}_{iq}\mathbf{x}-{\mathbf{b}}_{iq}\ge \mathbf{o}\end{array}$$
中的$\mathbf{H}$、$\mathbf{c}$、${\mathbf{A}}_{eq}$、${\mathbf{b}}_{eq}$、${\mathbf{A}}_{iq}$和${\mathbf{b}}_{iq}$，x1表示问题的初始可行解${\mathbf{x}}_1$。其中，Aeq，beq，Aiq和biq是命名参数，缺省值为无类型常量None。
函数体内，第4~{}18行将表示${\mathbf{A}}_{eq}$和${\mathbf{A}}_{iq}$的Aeq和Aiq叠加成$\mathbf{A}=\left(\begin{array}{c}{\mathbf{A}}_{eq}\\ {\mathbf{A}}_{iq}\end{array}\right)$，赋予A。将表示${\mathbf{b}}_{eq}$和${\mathbf{b}}_{iq}$的beq和biq叠加成$\mathbf{b}=\left(\begin{array}{c}{\mathbf{b}}_{eq}\\ {\mathbf{b}}_{iq}\end{array}\right)$，赋予b。记录下等式约束个数l和不等式约束个数m。
第19行定义问题(4)的目标函数的梯度$\nabla f(\mathbf{x})=\mathbf{Hx}+\mathbf{c}$赋予gk。第20、21行设置等式约束下标集E和不等式约束下标集I。第22~24初始化迭代次数k，迭代点xk，零向量b1。第25初始化xk处的有效集Ik，第26行初始化子问题的等式约束下标集Sk。
第28~51行的while循环，执行算法的迭代运算。其中，第29行调用博文《最优化方法Python计算：二次规划的拉格朗日算法》定义的拉格朗日算法函数qlag，求解
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{d}}^{\top }\mathbf{Hd}+\nabla f({\mathbf{x}}_k{)}^{\top }\mathbf{d}\\ \text{s.t. }{\mathbf{a}}_i\mathbf{d}=0,i\in I_k\end{array},$$
以求搜索方向${\mathbf{d}}_k$及其对应的拉格朗日乘子${\mathbf{\lambda }}_k$，赋予dk和lamdk。第30~41行处理当搜索方向${\mathbf{d}}_k\ne \mathbf{o}$时，按
α k = min ⁡ { 1 , α ^ k } \alpha_k=\min\{1,\hat{\alpha}_k\} αk​=min{1,α^k​}
其中，${\hat{\alpha }}_k=\min \left(\frac{b_i-{\mathbf{a}}_i{\mathbf{x}}_k}{{\mathbf{a}}_i{\mathbf{d}}_k}|i\in \overline{I_k},{\mathbf{a}}_i{\mathbf{d}}_k<0\right)$，计算搜索步长${\alpha }_k$，赋予alpha。以xk+alpha*dk更新迭代点xk，并更新有效集Ik和子问题等式约束下标集Sk。第42~51行处理所得搜索方向${\mathbf{d}}_k=\mathbf{o}$的情形。根据拉格朗日乘子${\mathbf{\lambda }}_k$是否存在负值元素，或更新有效集并进入下一轮迭代（第44~49行）或完成迭代置标志opt为True（第51行）。第52~54行返回最优解xk及最优目标近似值。
例1： 用qpact函数求二次规划
$$\{\begin{array}{l}\text{minimize}{x}_1^2+2x_2^2+x_3^2-2x_1{x}_2+x_3\\ \text{s.t. }{x}_1+x_2+x_3=4\\ 2x_1-x_2+x_3=2\end{array},$$
给定初始可行解${\mathbf{x}}_1=\left(\begin{array}{c}2\\ 2\\ 0\end{array}\right)$。
解：本题中，
$$\mathbf{H}=\left(\begin{array}{ccc}2& -2& 0\\ -2& 4& 0\\ 0& 0& 2\end{array}\right),{\mathbf{A}}_{eq}=\left(\begin{array}{ccc}1& 1& 1\\ 2& -1& 1\end{array}\right),{\mathbf{b}}_{eq}=\left(\begin{array}{c}4\\ 2\end{array}\right),\mathbf{c}=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$$
下列代码完成计算。

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.977272727272721
 nit: 1
   x: array([21/11, 43/22, 3/22])

仅迭代1次即算得最优解${\mathbf{x}}_0=\left(\begin{array}{c}\frac{21}{11}\\ \frac{43}{22}\\ \frac{3}{22}\end{array}\right)$。
例2：用qpact函数求解二次规划
$$\{\begin{array}{l}\text{minimize}{x}_1^2-x_1{x}_2+2x_2^2-x_1-10x_2\\ \text{s.t. }-3x_1-2x_2\ge -6\\ x_1,x_2\ge 0\end{array},$$
给定初始可行解${\mathbf{x}}_1=\mathbf{o}$。
解：本题中
$$\mathbf{H}=\left(\begin{array}{cc}2& -1\\ -1& 4\end{array}\right),\mathbf{c}=\left(\begin{array}{c}-1\\ -10\end{array}\right),{\mathbf{x}}_1=\left(\begin{array}{c}0\\ 0\end{array}\right),{\mathbf{A}}_{iq}=\left(\begin{array}{cc}-3& -2\\ 1& 0\\ 0& 1\end{array}\right),{\mathbf{b}}_{iq}=\left(\begin{array}{c}-6\\ 0\\ 0\end{array}\right).$$
下列代码利用这些数据完成本例计算。

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.75
 nit: 3
   x: array([1/2, 9/4])

这意味着，经过3此迭代，有效集算法算得本问题最优解${\mathbf{x}}_0=\left(\begin{array}{c}\frac{1}{2}\\ \frac{9}{4}\end{array}\right)$，最优值$f({\mathbf{x}}_0)\approx -13.75$。
例3：用qpact函数求解二次规划
$$\{\begin{array}{l}\text{minimize}{x}_1^2+x_1{x}_2+2x_2^2+x_3^2-6x_1-2x_2-12x_3\\ \text{s.t. }{x}_1+x_2+x_3-2=0\\ x_1-2x_2+3\ge 0\\ x_1,x_2,x_3\ge 0\end{array},$$
给定初始可行解${\mathbf{x}}_1=\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right)$。
{\heiti{解}}：本问题中的数据矩阵表示为
$$\begin{array}{l}\mathbf{H}=\left(\begin{array}{ccc}2& 1& 0\\ 1& 4& 0\\ 0& 0& 2\end{array}\right),\mathbf{c}=\left(\begin{array}{c}-6\\ -2\\ -12\end{array}\right)\\ {\mathbf{A}}_{eq}=\left(\begin{array}{ccc}1& 1& 1\end{array}\right),{\mathbf{b}}_{eq}=(2)\\ {\mathbf{A}}_{iq}=\left(\begin{array}{ccc}1& -2& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),{\mathbf{b}}_{iq}=\left(\begin{array}{c}-3\\ 0\\ 0\\ 0\end{array}\right)\end{array}$$
下列代码完成计算

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.0
 nit: 2
   x: array([0, 0, 2])

即2次迭代，算得最优解为${\mathbf{x}}_0=\left(\begin{array}{c}0\\ 0\\ 2\end{array}\right)$，最优值为$f({\mathbf{x}}_0)=-20$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！