控制理论：MindOpt Python API求解模型预测控制问题（Model Predictive Control）

1. MPC应用背景

模型预测控制（Model Predictive Control，MPC），也被称为滚动时域控制（Receding Horizon Control），是优化和控制两个领域的交叉。其基本思想是在每个需要控制的时刻预测受控系统在未来有限时间范围内行为，并基于预测结果计算当前时刻的最佳控制输入，在确保满足系统约束的同时，最小化成本函数。这个过程会反复进行，是一个滚动的动态过程。

以开车举例，假设司机行驶在路上，发现前面车辆踩刹车，车尾亮红灯。那在此刻有限时间内对系统未来行为的预测是前车减速，需要在满足自身车辆加速度约束的情况下，计算踩刹车的力度。此时，若踩刹车力度过大踩的太急容易发生后车追尾事故，因此需要控制踩刹车的力度，该力度可被视为最佳控制输入。在不变道的情况下，该过程是反复出现的，需要根据每一次道路上的情况对系统行为进行预测，从而动态的调整踩刹车/油门力度（油门可视为负的刹车）。

综上，MPC控制本质上可看成，通过构建模型预测未来有限时间内的变化，在满足约束的情况下，动态决策当前系统的最优控制量。由于预测是针对未来有限时间范围进行的，在每一个控制时刻都需要针对当前时刻的状态去迭代求解。

---

2. 问题描述

给定一个线性时不变动态系统（Linear Time-invariant Dynamical System），考虑通过模型预测控制实现对该系统的最优控制。假设该系统的控制时刻集合为$T$，系统在$t$时刻状态由变量$x_t\in R^{n_x}$描述，系统输入由变量$u_t\in R^{n_u}$描述。通过推导可将该问题转化为一个二次规划模型，目标代价函数仅跟初始条件和输入量有关。在优化模型创建章节将对二次规划模型进行详细阐述。

3. 优化模型构建

下面我们构建二次规划模型求解MPC控制问题，分别定义该数学模型的决策变量、目标函数以及约束等。

具体的集合、参数与决策变量定义如下。

符号	含义
$T$	控制时刻集合，任意控制时刻用$indext\in T$表示；
$Q$	目标函数系数矩阵，状态代价矩阵，对称矩阵；
$R$	目标函数系数矩阵，控制代价矩阵，对称矩阵；
$A$	系数矩阵，状态转移系数，特征值小于1；
$B$	系数矩阵，状态转移系数；
$x_{min}$	参数，状态取值下限；
$x_{max}$	参数，状态取值上限；
$u_{min}$	参数，输入取值下限；
$u_{max}$	参数，输入取值上限；
$x_{cur}$	参数，初始状态取值；
$n_x$	参数，状态变量的维度；
$n_u$	参数，输入变量的维度；
决策变量	含义
$x_t$	连续变量，表示$t$时刻状态取值；
$u_t$	连续变量，表示$t$时刻输入取值；

• 目标函数：代价函数；

$minz={\sum }_{t=0}^{T-1}{x}_t^{T}Q{x}_t+u_t^{T}R{u}_t+x_T^{T}Q{x}_T$

• 约束1：状态转移约束，$t+1$时刻的状态变量由$t$时刻状态变量与$t$时刻的系统输入变量共同决定；

x t + 1 = A x t + B u t ∀ t ∈ T − { T } x_{t+1} = Ax_t + Bu_t \quad \forall t \in T-\{T\} xt+1​=Axt​+But​∀t∈T−{T}

• 约束2：状态取值上下限约束；

$x_{min}\le x_t\le x_{max}\forall t\in T$

• 约束3：输入取值上下限约束；

$u_{min}\le u_t\le u_{max}\forall t\in T$

• 约束4：初始状态取值约束；

$x_0=x_{cur}$

---

在后续内容中，分别展示对上述MPC问题数学模型进行单次建模求解，以及迭代建模求解的代码。

其中，章节4中展示了单次建模代码，章节5展示了迭代求解建模代码。

4. MindOpt Python API模型输入和求解

通过使用MindOpt API，可以通过不同类型的编程语言对数学规划模型进行建模求解，比如C、Python、Java等，具体的操作见用户手册MindOpt Python API链接。以下我们通过使用MindOpt Python API，以按行输入的形式来对上述QP问题进行建模与求解。

首先，引入建模用到的Python包：

import numpy as np
import scipy as sp
from scipy import sparse
import pandas as pd
from mindoptpy import *

OpenBLAS WARNING - could not determine the L2 cache size on this system, assuming 256k
OpenBLAS WARNING - could not determine the L2 cache size on this system, assuming 256k

之后，输入建模中运用到的参数，并对参数矩阵进行预处理。

其中，矩阵Q是目标函数中的状态变量系数对称矩阵，矩阵R是目标函数中的输入变量系数对称矩阵。矩阵A与矩阵B是系统模型系数矩阵。除此之外，输入的参数包括初始状态参数向量，状态变量上下限参数，输入变量上下限参数等。

# Coefficient matrix
A = sparse.csc_matrix([
 [ 9.38140422e-01, -7.60325726e-03, -8.46189503e-02,  1.00709617e-01,
  -1.10581661e-01, -5.48517782e-02, -1.19182527e-01,  9.20903867e-02,
  -8.13588333e-02, -5.32070948e-02,  1.29533953e-01, -9.87224752e-02],
 [ 7.16523283e-02,  9.28612532e-01,  2.42877761e-02, -3.59338778e-02,
  -9.44847432e-02,  7.25071608e-02, -1.30168031e-01,  4.68817358e-03,
   3.14799623e-02, -8.82754918e-02, -1.28610493e-02,  2.40896182e-02],
 [ 6.07909138e-02, -3.79483900e-02,  9.19755623e-01,  1.70227210e-02,
  -4.08411644e-02, -5.81220168e-02,  1.27332335e-02,  6.35131066e-02,
  -7.45012238e-02, -5.01856212e-03, -5.27963356e-02, -7.71163246e-02],
 [ 5.13533756e-02, -4.51163548e-02,  1.65354520e-01,  9.01457399e-01,
  -3.83313492e-02, -2.14163200e-03,  1.16619058e-01,  3.94289072e-03,
   7.14068192e-02,  8.15587042e-02,  1.87321777e-02, -8.58947531e-02],
 [-4.33623776e-03,  7.27227727e-02, -6.95740426e-03, -2.51122321e-02,
   1.01251572e+00,  7.67482873e-02,  3.62972504e-03,  1.14530753e-02,
   4.29389660e-02, -8.68872513e-03, -3.42707613e-02,  1.11016572e-02],
 [-2.14918207e-02,  8.10484863e-03, -6.57621215e-02,  2.25303519e-02,
   3.96897204e-02,  8.14805249e-01,  8.40466742e-02,  2.62653448e-02,
  -1.48686414e-01,  7.57462329e-02,  1.06799488e-02, -9.15667238e-02],
 [-8.64200350e-02,  7.27862778e-02, -5.59151529e-02, -2.05061517e-02,
   1.89439169e-01, -5.01285381e-03,  9.61857585e-01,  1.44976402e-02,
  -1.01510017e-01,  7.61703988e-03, -1.12474461e-02,  1.16469238e-01],
 [-6.93556115e-04, -7.72577198e-03,  5.68327080e-03, -9.76048766e-02,
   1.05252420e-01,  5.69580141e-02, -6.35334019e-02,  8.78118282e-01,
  -7.14207427e-02,  5.07488700e-02, -6.69718832e-02,  5.24863217e-02],
 [ 4.95802757e-02, -5.64816811e-02, -7.99671506e-02,  9.76317500e-02,
  -3.99783452e-02, -1.53551987e-01,  7.86414044e-02,  5.04390435e-02,
   8.44393950e-01,  2.57175522e-02,  2.65963742e-03, -5.06924976e-02],
 [ 1.09353418e-02,  8.22376384e-02, -6.62280502e-03, -2.74759885e-02,
  -5.97409376e-02,  6.65583342e-02, -1.49478036e-02,  1.74621538e-02,
   8.20676675e-02,  9.56606140e-01, -2.00780101e-02, -4.50562645e-02],
 [-1.37047251e-02,  5.77684191e-02,  1.82589599e-02, -9.39575095e-02,
   1.33644617e-01,  2.41672011e-02, -2.57816950e-02,  8.47327486e-02,
  -2.44082896e-03, -4.92168263e-02,  9.39514778e-01,  1.01267032e-01],
 [ 3.67429627e-02, -5.47835358e-02, -6.75801432e-02,  6.65385440e-02,
  -6.97549932e-02,  3.54784024e-02, -1.79539375e-01, -1.18970587e-01,
  -1.26535985e-01, -1.08899804e-02, -1.28202973e-02,  9.95029868e-01]
])
B = sparse.csc_matrix([
 [-2.24360505,  1.23346365, -0.08718269, -0.65867254, -0.59175896,  1.10527788],
 [-1.0237667, -1.06206018,  1.36152413, -0.21642054,  0.41134999, 0.99929514],
 [ 0.40566331, -0.94471656, -0.18820403,  0.18618359, -1.06582946, -1.68116081],
 [-1.21377262, -1.22278943,  0.05721832, -0.05359124,  1.665432,    0.8250249 ],
 [ 0.69572001,  1.48726486,  1.94257856,  1.03655593, -0.70578538,  0.03685261],
 [ 0.79515426,  1.39900998, -1.44463614, -1.4922319,  -0.07020333,  1.04563105],
 [-0.52454659, -2.12389802,  0.87192919,  0.50470015,  1.1183585,   1.3051162 ],
 [-0.14954945,  1.45139811,  0.26399847,  0.61619104, -0.10564964, -0.04489195],
 [-1.21396425, -2.14118422,  0.23671791,  1.13060713,  0.29053468,  0.14718403],
 [-0.73695489,  0.01549698, -0.09137905,  0.68290566, -0.0933051,   0.08755412],
 [-0.92919747, -0.82959582,  2.40329256,  0.21520008, -0.88681758, -0.98510142],
 [ 1.07819139, -1.78632351, -0.57792416,  0.64875391, -0.72312702, -0.18435067]
])

nx, nu = 12, 6

# Upper bound & lower bound of variables
umin = np.array([-0.10621128,-0.2806848, -0.77621463, -0.8901663,  -0.51059931, -0.71570793])
umax = -umin
xmin = np.array([-1.07777887, -1.61530423, -1.83780387, -1.6960054,  -1.55505936, -1.86118965,
 -1.74532151, -1.26454931, -1.99976912, -1.5737319,  -1.63905897, -1.4896996 ])
xmax = -xmin

# Objective function 
Q = np.array([[ 1.57078007e+00, -5.78298110e-01,  1.93921610e-01,  4.65120772e-01,
   4.94597667e-01, -4.98092806e-01,  2.19148227e-01, -6.46375537e-01,
  -1.48233557e-01, -3.66068757e-01, -1.90480898e-01, -1.66313639e-02],
 [-5.78298110e-01,  3.96159732e+00,  1.12371922e-01, -1.01921812e+00,
  -1.81435516e+00, -2.21013948e-01, -3.00613156e+00,  1.05914941e+00,
   7.34762833e-01, -3.42910086e-01, -1.69944688e+00, -3.53296139e-02],
 [ 1.93921610e-01,  1.12371922e-01,  1.28074366e-01,  9.35286575e-02,
   7.57506315e-02, -1.45982736e-01, -5.06976522e-02, -5.10636822e-02,
  -5.38402084e-02, -1.00903209e-01, -2.44899837e-01,  7.83643369e-04],
 [ 4.65120772e-01, -1.01921812e+00,  9.35286575e-02,  5.44685736e-01,
   7.98114795e-01, -2.15324795e-01,  5.66797260e-01, -8.56484512e-01,
  -9.04553190e-02, -7.54041779e-02,  1.76940708e-01, -1.16593157e-01],
 [ 4.94597667e-01, -1.81435516e+00,  7.57506315e-02,  7.98114795e-01,
   2.29184674e+00, -1.16916496e-01,  1.27818732e+00, -1.33451641e+00,
  -9.97041479e-02,  6.95565274e-02,  5.06620933e-01, -2.84192560e-01],
 [-4.98092806e-01, -2.21013948e-01, -1.45982736e-01, -2.15324795e-01,
  -1.16916496e-01,  1.33133012e+00,  5.78028958e-01,  6.67124759e-01,
  -1.89241194e-01,  5.45216087e-01,  6.75826169e-01,  1.27294772e-01],
 [ 2.19148227e-01, -3.00613156e+00, -5.06976522e-02,  5.66797260e-01,
   1.27818732e+00,  5.78028958e-01,  4.10574361e+00, -4.42980591e-02,
  -9.97013988e-01,  5.85841617e-01,  1.69327526e+00,  2.13040400e-01],
 [-6.46375537e-01,  1.05914941e+00, -5.10636822e-02, -8.56484512e-01,
  -1.33451641e+00,  6.67124759e-01, -4.42980591e-02,  2.96841592e+00,
  -6.38029875e-01,  2.31628466e-01,  2.98128801e-01,  7.27663907e-01],
 [-1.48233557e-01,  7.34762833e-01, -5.38402084e-02, -9.04553190e-02,
  -9.97041479e-02, -1.89241194e-01, -9.97013988e-01, -6.38029875e-01,
   5.78359880e-01, -1.19071026e-01, -5.25589347e-01, -3.28725354e-01],
 [-3.66068757e-01, -3.42910086e-01, -1.00903209e-01, -7.54041779e-02,
   6.95565274e-02,  5.45216087e-01,  5.85841617e-01,  2.31628466e-01,
  -1.19071026e-01,  4.12306092e-01,  5.43418830e-01,  2.83159994e-02],
 [-1.90480898e-01, -1.69944688e+00, -2.44899837e-01,  1.76940708e-01,
   5.06620933e-01,  6.75826169e-01,  1.69327526e+00,  2.98128801e-01,
  -5.25589347e-01,  5.43418830e-01,  1.84403232e+00,  2.33772964e-01],
 [-1.66313639e-02, -3.53296139e-02,  7.83643369e-04, -1.16593157e-01,
  -2.84192560e-01,  1.27294772e-01,  2.13040400e-01,  7.27663907e-01,
  -3.28725354e-01,  2.83159994e-02,  2.33772964e-01,  2.84744003e-01]]) * 0.5
R = 0.1 * np.eye(nu) * 0.5

Q1, R1 = Q.tolist(), R.tolist()

QN = Q1


# Initial states
x0 = np.array([0.41965988, 0.60494802, 0.63916977, 0.25054046, 0.1330799, -0.40651584,
  0.85809431, 0.37893636, 0.00327209, 0.42855133, -0.61582717, 0.62719014])

# Prediction horizon
T = 10

通过MindOpt Python API创建优化模型，并在优化模型中逐步添加变量、约束、目标函数等，在完成添加后通过MindOpt求解器对优化模型进行求解。

• 通过调用model.addVar()来添加模型中涉及的优化变量，将约束2与约束3转化为变量$x_t$与$u_t$的上下界，减少了约束的添加。
• 通过调用model.addConstr()来添加模型中涉及的线性等式约束1&约束4；其中，通过quicksum()能够以向量形式添加约束，加快求解速度。
• 为模型设置目标函数时，首先使用类QuadExpr()创建一个二次表达式， 利用QuadExpr中的方法QuadExpr.addTerms()输入目标函数的二次部分，最后再用Model.setObjective()来设置目标函数并将问题设置为最小化。
• 数学模型建立完成后，再调用Model.optimize()求解优化问题。

# Create model
model = Model("MPC_01")

# Add variables.
x = dict()
u = dict()

# Add variables x
for t in range(T + 1):
    for i in range(nx):
        x[t, i] = model.addVar(xmin[i], xmax[i], 0.0, 'C', "x_" + str(t) + "_" + str(i))

# Add variables u
for t in range(T):
    for i in range(nu):
        u[t, i] = model.addVar(umin[i], umax[i], 0.0, 'C', "u_" + str(t) + "_" + str(i))
        
# Add linear constraints
# constr 4:
for i in range(nx):
    model.addConstr(x[0, i] == (x0[i], x0[i]), name = "constr4_" + str(i))
    
# constr 1:
for t in range(T):
    for i in range(nx):
        model.addConstr(x[t+1, i] == quicksum(A[i, jx] * x[t, jx] for jx in range(nx)) + 
                        quicksum(B[i, ju] * u[t, ju] for ju in range(nu)), name = "constr1_" + str(t) + "_" + str(i))
    
# Add Objective Function
obj = QuadExpr()

for t in range(T):
    for i in range(nx):
        obj.addTerms(Q1[i], [x[t, i] for j in range(nx)], [x[t, j] for j in range(nx)])
    
    for i in range(nu):
        obj.addTerms(R1[i], [u[t, i] for j in range(nu)], [u[t, j] for j in range(nu)])

for i in range(nx):
    obj.addTerms(QN[i], [x[T, i] for j in range(nx)], [x[T, j] for j in range(nx)])

model.setObjective(obj, MDO.MINIMIZE)

# Solve Model
model.optimize()

# Print optimal objective function
if model.status != MDO.OPTIMAL:
    raise ValueError('MindOpt did not solve the problem!')

print('-' * 50)
print(f"Optimal objective value is: {model.objval}")

运行结果如下：

Start license validation (current time : 30-JUL-2024 16:26:20).
License validation terminated. Time : 0.004s

Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [3.3e-03,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.002s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.59415202e+01 -2.51628148e+02 5.3e-01 3.4e-01 8.0e+00 9.5e-02   0.02s
    1 +3.65894235e+00 -7.57131305e+00 1.5e-03 9.5e-04 3.1e+00 7.4e-03   0.03s
    2 +3.59251638e+00 +3.34434913e+00 4.0e-06 2.0e-04 8.4e-02 1.9e-04   0.03s
    3 +3.58895136e+00 +3.58686168e+00 1.1e-08 4.5e-06 6.5e-04 1.6e-06   0.03s
    4 +3.58890782e+00 +3.58890276e+00 2.0e-11 1.8e-07 1.4e-06 1.6e-09   0.03s
    5 +3.58890780e+00 +3.58890779e+00 5.6e-14 4.9e-10 3.9e-09 4.3e-12   0.03s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 3.5889078011199E+00
 - Dual objective     : 3.5889077872126E+00
 - Num. threads       : 2
 - Num. iterations    : 5
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.016s

--------------------------------------------------
Optimal objective value is: 3.5889078011198436

---

5.迭代求解MPC问题

通过迭代求解该QP模型可实现通过求解MPC问题对动态系统的控制,具体方式为：

在每次完成QP问题的求解后可得到当前阶段输入变量$u_0$的取值，通过以下公式修正当前状态$x_{cur}$。

$x_{cur}=A{x}_0+B{u}_0$

不断迭代求解QP问题，实现对动态系统的持续控制，具体流程如下。

# Create model
model = Model("MPC_01")

# Add variables.
x = dict()
u = dict()

# Add variables x
for t in range(T + 1):
    for i in range(nx):
        x[t, i] = model.addVar(xmin[i], xmax[i], 0.0, 'C', "x_" + str(t) + "_" + str(i))

# Add variables u
for t in range(T):
    for i in range(nu):
        u[t, i] = model.addVar(umin[i], umax[i], 0.0, 'C', "u_" + str(t) + "_" + str(i))
        
# Add linear constraints
# constr 4:
for i in range(nx):
    model.addConstr(x[0, i] == (x0[i], x0[i]), name = "constr4_" + str(i))
    
# constr 1:
for t in range(T):
    for i in range(nx):
        model.addConstr(x[t+1, i] == quicksum(A[i, jx] * x[t, jx] for jx in range(nx)) + 
                        quicksum(B[i, ju] * u[t, ju] for ju in range(nu)), name = "constr1_" + str(t) + "_" + str(i))
    
# Add Objective Function
obj = QuadExpr()

for t in range(T):
    for i in range(nx):
        obj.addTerms(Q1[i], [x[t, i] for j in range(nx)], [x[t, j] for j in range(nx)])
    
    for i in range(nu):
        obj.addTerms(R1[i], [u[t, i] for j in range(nu)], [u[t, j] for j in range(nu)])

for i in range(nx):
    obj.addTerms(QN[i], [x[T, i] for j in range(nx)], [x[T, j] for j in range(nx)])

model.setObjective(obj, MDO.MINIMIZE)

result = list()

# Number of simulation
Nsim = 5

for _ in range(Nsim):
    # Solve model
    model.optimize()

    # Print optimal objective function
    if model.status != MDO.OPTIMAL:
        raise ValueError('MindOpt did not solve the problem!')

    print(f"Optimal objective value is: {model.objval}")
    
    # store obj
    result.append(model.objVal)
    
    # update initial states
    u0 = np.array([u[0, uidx].x for uidx in range(nu)])

    x0 = A@x0 + B@u0

    
    # Create model
    model = Model("MPC_" + str(_))

    # Add variables.
    x = dict()
    u = dict()

    # Add variables x
    for t in range(T + 1):
        for i in range(nx):
            x[t, i] = model.addVar(xmin[i], xmax[i], 0.0, 'C', "x_" + str(t) + "_" + str(i))

    # Add variables u
    for t in range(T):
        for i in range(nu):
            u[t, i] = model.addVar(umin[i], umax[i], 0.0, 'C', "u_" + str(t) + "_" + str(i))

    # Add linear constraints
    # constr 4:
    for i in range(nx):
        model.addConstr(x[0, i] == (x0[i], x0[i]), name = "constr4_" + str(i))

    # constr 1:
    for t in range(T):
        for i in range(nx):
            model.addConstr(x[t+1, i] == quicksum(A[i, jx] * x[t, jx] for jx in range(nx)) + 
                            quicksum(B[i, ju] * u[t, ju] for ju in range(nu)), name = "constr1_" + str(t) + "_" + str(i))

    # Add Objective Function
    obj = QuadExpr()

    for t in range(T):
        for i in range(nx):
            obj.addTerms(Q1[i], [x[t, i] for j in range(nx)], [x[t, j] for j in range(nx)])

        for i in range(nu):
            obj.addTerms(R1[i], [u[t, i] for j in range(nu)], [u[t, j] for j in range(nu)])

    for i in range(nx):
        obj.addTerms(QN[i], [x[T, i] for j in range(nx)], [x[T, j] for j in range(nx)])

    model.setObjective(obj, MDO.MINIMIZE)

运行结果如下：

Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [3.3e-03,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.59415202e+01 -2.51628148e+02 5.3e-01 3.4e-01 8.0e+00 9.5e-02   0.01s
    1 +3.65894235e+00 -7.57131305e+00 1.5e-03 9.5e-04 3.1e+00 7.4e-03   0.01s
    2 +3.59251638e+00 +3.34434913e+00 4.0e-06 2.0e-04 8.4e-02 1.9e-04   0.01s
    3 +3.58895136e+00 +3.58686168e+00 1.1e-08 4.5e-06 6.5e-04 1.6e-06   0.01s
    4 +3.58890782e+00 +3.58890276e+00 2.0e-11 1.8e-07 1.4e-06 1.6e-09   0.02s
    5 +3.58890780e+00 +3.58890779e+00 5.6e-14 4.9e-10 3.9e-09 4.3e-12   0.02s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 3.5889078011198E+00
 - Dual objective     : 3.5889077872126E+00
 - Num. threads       : 2
 - Num. iterations    : 5
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.007s

Optimal objective value is: 3.588907801119851
Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [1.3e-02,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.34176239e+01 -2.56954966e+02 4.9e-01 3.4e-01 8.7e+00 1.0e-01   0.01s
    1 +2.31477425e+00 -9.54340882e+00 1.4e-03 9.4e-04 5.1e+00 7.8e-03   0.01s
    2 +2.24708454e+00 +1.96341892e+00 3.7e-06 2.6e-05 1.5e-01 1.9e-04   0.01s
    3 +2.24591415e+00 +2.24248149e+00 1.0e-08 1.5e-05 1.8e-03 2.8e-06   0.01s
    4 +2.24591337e+00 +2.24590366e+00 2.8e-11 4.0e-08 5.1e-06 7.8e-09   0.01s
    5 +2.24591337e+00 +2.24591334e+00 7.7e-14 1.1e-10 1.4e-08 2.1e-11   0.01s
    6 +2.24591337e+00 +2.24591337e+00 2.3e-16 3.0e-13 3.9e-11 5.9e-14   0.01s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 2.2459133670016E+00
 - Dual objective     : 2.2459133669277E+00
 - Num. threads       : 2
 - Num. iterations    : 6
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.007s

Optimal objective value is: 2.2459133670015845
Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [1.2e-02,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.27679533e+01 -2.53526951e+02 4.8e-01 3.4e-01 8.7e+00 9.7e-02   0.01s
    1 +1.91978324e+00 -9.55624291e+00 1.3e-03 9.4e-04 6.0e+00 7.5e-03   0.01s
    2 +1.84891367e+00 +1.67169445e+00 3.6e-06 5.0e-05 1.1e-01 1.2e-04   0.01s
    3 +1.84780269e+00 +1.84587415e+00 9.9e-09 6.2e-06 1.2e-03 1.4e-06   0.01s
    4 +1.84780217e+00 +1.84779678e+00 2.7e-11 1.7e-08 3.2e-06 4.0e-09   0.01s
    5 +1.84780217e+00 +1.84780215e+00 7.5e-14 4.7e-11 8.9e-09 1.1e-11   0.01s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 1.8478021683519E+00
 - Dual objective     : 1.8478021534841E+00
 - Num. threads       : 2
 - Num. iterations    : 5
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.006s

Optimal objective value is: 1.8478021683519037
Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [5.1e-03,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.24282214e+01 -2.49583321e+02 4.8e-01 3.4e-01 8.7e+00 9.3e-02   0.01s
    1 +1.69129612e+00 -9.32067856e+00 1.3e-03 9.4e-04 6.5e+00 7.2e-03   0.01s
    2 +1.62970160e+00 +1.47852487e+00 3.6e-06 4.9e-05 1.1e-01 1.1e-04   0.01s
    3 +1.62891312e+00 +1.62771692e+00 1.0e-08 4.1e-06 8.2e-04 9.0e-07   0.01s
    4 +1.62891296e+00 +1.62890965e+00 2.7e-11 1.1e-08 2.3e-06 2.5e-09   0.02s
    5 +1.62891296e+00 +1.62891295e+00 7.5e-14 3.1e-11 6.2e-09 6.9e-12   0.02s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 1.6289129626863E+00
 - Dual objective     : 1.6289129535508E+00
 - Num. threads       : 2
 - Num. iterations    : 5
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.008s

Optimal objective value is: 1.6289129626861671
Model summary.
 - Num. variables     : 192
 - Num. constraints   : 132
 - Num. nonzeros      : 2292
 - Bound range        : [2.6e-04,2.0e+00]
 - Quad. bound range  : [1.1e-01,2.0e+00]
 - Objective range    : [0.0e+00,0.0e+00]
 - Quad. obj. range   : [7.8e-04,4.1e+00]
 - Matrix range       : [6.9e-04,2.4e+00]

Presolver started.
Presolver terminated. Time : 0.001s

Interior point method started.
 Iter         PrimObj         DualObj PrimFea DualFea  GapFea      Mu   Time
    0 +3.20626216e+01 -2.54853451e+02 4.9e-01 3.4e-01 8.9e+00 9.8e-02   0.01s
    1 +1.50527817e+00 -1.01054960e+01 1.3e-03 9.4e-04 7.7e+00 7.6e-03   0.01s
    2 +1.45014309e+00 +1.29394541e+00 3.7e-06 4.7e-05 1.3e-01 1.1e-04   0.01s
    3 +1.44946501e+00 +1.44844590e+00 1.0e-08 1.9e-06 7.5e-04 7.4e-07   0.01s
    4 +1.44946493e+00 +1.44946211e+00 2.8e-11 5.3e-09 2.1e-06 2.0e-09   0.01s
    5 +1.44946493e+00 +1.44946492e+00 7.7e-14 1.5e-11 5.7e-09 5.6e-12   0.01s
Terminated.
 - Method             : Interior point method.
 - Primal objective   : 1.4494649256440E+00
 - Dual objective     : 1.4494649178793E+00
 - Num. threads       : 2
 - Num. iterations    : 5
 - Solver details     : Solver terminated with a primal/dual optimal status.

Interior point method terminated. Time : 0.006s

Optimal objective value is: 1.4494649256439367

---

最后，将存储的求解结果通过pandas输出为csv文件：

# column name
columns = ['Time Index', 'Objective Function']

# create dataframe
df_Output = pd.DataFrame(columns=columns)
for t in range(len(result)):
    df_Output.loc[len(df_Output)] = [t, result[t]]

# output
df_Output.to_csv('Result.csv', index=False)

6.运行结果

迭代重复计算QP模型，并记录计算结果。记录的解被存储在[Result.csv]文件中，具体结果如下所示。

Time Index	Objective Function
0	3.5889078011198485
1	2.2459133670015876
2	1.8478021683502268
3	1.6289129620867264
4	1.4494649240215713