文章目录
- 一、问题描述
- 二、数学推导
-
- 1. 目标函数处理
- 2. 约束条件处理
- 三、代码编写
一、问题描述
已知:
m
i
n
(
x
1
−
1
)
2
+
(
x
2
−
2
)
2
s
.
t
.
0
⩽
x
1
⩽
1.5
,
1
⩽
x
2
⩽
2.5
min(x_1-1)^2+(x_2-2)^2 \\qquad s.t. \\ \\ 0 \\leqslant x_1 \\leqslant 1.5,\\ \\ 1 \\leqslant x_2 \\leqslant 2.5
min(x1−1)2+(x2−2)2s.t. 0⩽x1⩽1.5, 1⩽x2⩽2.5 目标函数为二元二次函数,可行域为线性、凸集,此为二次规划问题,可将其转换成二次规划表达式再进行求解。相关数学概念参考另一篇: 最优化问题基础理论概述。
二、数学推导
1. 目标函数处理
f
(
x
1
,
x
2
)
=
(
x
1
−
1
)
2
+
(
x
2
−
2
)
2
=
x
1
2
+
x
2
2
−
2
x
1
−
4
x
2
+
C
f(x_1, x_2)=(x_1-1)^2+(x_2-2)^2 =x_1^2+x_2^2-2x_1-4x_2+C
f(x1,x2)=(x1−1)2+(x2−2)2=x12+x22−2x1−4x2+C
其中,常数项用
C
C
C表示;
令,
X
=
[
x
1
x
2
]
X=\\left[ \\begin{matrix} x_1 \\\\ x_2 \\end{matrix} \\right]
X=[x1x2],则
f
(
x
1
,
x
2
)
=
[
x
1
,
x
2
]
[
x
1
,
x
2
]
T
+
[
−
2
,
−
4
]
[
x
1
,
x
2
]
T
=
X
T
X
+
[
−
2
,
−
4
]
X
=
1
2
X
T
[
2
0
0
2
]
X
+
[
−
2
−
4
]
T
X
=
1
2
X
T
P
X
+
Q
T
X
\\begin{aligned} f(x_1, x_2) &= [x_1, x_2][x_1, x_2]^T+[-2, -4][x_1, x_2]^T \\\\[2ex] &= X^TX+[-2, -4]X \\\\[2ex] &=\\frac{1}{2} X^T \\left[\\begin{matrix} 2 &0 \\\\ 0&2 \\end{matrix} \\right] X+\\left[\\begin{matrix} -2 \\\\ -4 \\end{matrix} \\right]^TX \\\\[2ex] &= \\frac{1}{2} X^TPX+Q^TX \\end{aligned}
f(x1,x2)=[x1,x2][x1,x2]T+[−2,−4][x1,x2]T=XTX+[−2,−4]X=21XT[2002]X+[−2−4]TX=21XTPX+QTX
其中,
P
=
[
2
0
0
2
]
,
Q
=
[
−
2
−
4
]
P=\\left[\\begin{matrix} 2 &0 \\\\ 0&2 \\end{matrix} \\right],\\ Q=\\left[\\begin{matrix} -2 \\\\ -4 \\end{matrix} \\right]
P=[2002], Q=[−2−4] 关于为什么要写成
1
2
X
T
P
X
\\frac{1}{2} X^TPX
21XTPX 形式,因为此时
P
P
P 为目标函数的海塞矩阵,具体参看 此链接。
2. 约束条件处理
{
0
⩽
x
1
⩽
1.5
1
⩽
x
2
⩽
2.5
⟺
[
0
1
]
⩽
[
x
1
x
2
]
⩽
[
1.5
2.5
]
⟺
[
0
1
]
⩽
[
1
0
0
1
]
X
⩽
[
1.5
2.5
]
\\begin{aligned} \\left\\{ \\begin{array}{} 0 \\leqslant x_1 \\leqslant 1.5 \\\\ 1 \\leqslant x_2 \\leqslant 2.5 \\\\ \\end{array} \\right . \\quad \\Longleftrightarrow \\quad \\left[\\begin{matrix}0 \\\\1\\end{matrix} \\right] \\leqslant \\left[\\begin{matrix}x_1 \\\\x_2\\end{matrix} \\right] \\leqslant \\left[\\begin{matrix}1.5 \\\\2.5\\end{matrix} \\right] \\quad \\Longleftrightarrow \\quad \\left[\\begin{matrix}0 \\\\1\\end{matrix} \\right] \\leqslant \\left[\\begin{matrix}1&0 \\\\0&1\\end{matrix} \\right]X \\leqslant \\left[\\begin{matrix}1.5 \\\\2.5\\end{matrix} \\right] \\end{aligned}
{0⩽x1⩽1.51⩽x2⩽2.5⟺[01]⩽[x1x2]⩽[1.52.5]⟺[01]⩽[1001]X⩽[1.52.5] 令
L
B
=
[
0
1
]
,
A
=
[
1
0
0
1
]
,
U
B
=
[
1.5
2.5
]
L_B=\\left[\\begin{matrix}0 \\\\1\\end{matrix} \\right],\\ A=\\left[\\begin{matrix}1&0 \\\\0&1\\end{matrix} \\right],\\ U_B=\\left[\\begin{matrix}1.5 \\\\2.5\\end{matrix} \\right]
LB=[01], A=[1001], UB=[1.52.5] ,整理得约束条件如下:
L
B
⩽
A
X
⩽
U
B
L_B \\leqslant AX \\leqslant U_B
LB⩽AX⩽UB
三、代码编写
由步骤 二、数学推导 得到5个矩阵:
-
P
P
-
Q
Q
-
U
B
U_B
-
L
B
L_B
-
A
A
现在根据这5个矩阵进行代码编写,是使用osqp进行二次型规划问题构建及求解。
代码如下:
Eigen::SparseMatrix<double> P(2, 2); // P, 二次型矩阵
Eigen::VectorXd Q(2); // Q, 一次项向量
Eigen::SparseMatrix<double> A(2, 2); // 单位阵
Eigen::VectorXd lowerBound(2); // 下边界向量
Eigen::VectorXd upperBound(2); // 上边界向量
P.insert(0, 0) = 2.0;
P.insert(1, 1) = 2.0;
std::cout << "\\033[34m" << "P:" << std::endl
<< P << "\\033[0m" << std::endl;
A.insert(0, 0) = 1.0;
A.insert(1, 1) = 1.0;
std::cout << "\\033[34m" << "A:" << std::endl
<< A << "\\033[0m" << std::endl;
Q << –2, –4;
std::cout << "\\033[34m" << "Q:" << std::endl
<< Q << "\\033[0m" << std::endl;
lowerBound << 0.0, 1.0;
upperBound << 1.5, 2.5;
// Step 1: 创建求解器
OsqpEigen::Solver solver;
// Step 2: 设置(提升求解速度)
solver.settings()->setVerbosity(false);
solver.settings()->setWarmStart(true);
// Step 3: 初始化(7部分)
solver.data()->setNumberOfVariables(2); // 变量数
solver.data()->setNumberOfConstraints(2); // 约束数
if (!solver.data()->setHessianMatrix(P)) // 海塞矩阵
{
return;
}
if (!solver.data()->setGradient(Q)) // Q矩阵
{
return;
}
if (!solver.data()->setLinearConstraintsMatrix(A)) // 线性约束矩阵A
{
return;
}
if (!solver.data()->setLowerBound(lowerBound)) // 下边界矩阵
{
return;
}
if (!solver.data()->setUpperBound(upperBound)) // 上边界矩阵
{
return;
}
if (!solver.initSolver())
{
return;
}
// Step 4:求解
Eigen::VectorXd QPSolution;
if (solver.solveProblem() != OsqpEigen::ErrorExitFlag::NoError)
{
return;
}
QPSolution = solver.getSolution();
std::cout << "\\033[1;32m" << "QPSolution:" << std::endl
<< QPSolution << "\\033[0m" << std::endl;
运行结果如下: 可见,当
x
1
=
1
,
x
2
=
2
x_1=1,\\ x_2=2
x1=1, x2=2 时目标函数取得最小。
评论前必须登录!
注册