使用osqp求解簡單二次規劃問題

文章目錄

一、問題描述
二、數學推導
- 1. 目標函數處理
- 2. 約束條件處理
三、代碼編寫

一、問題描述

已知：
$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$
目標函數為二元二次函數，可行域為線性、凸集，此為二次規劃問題，可將其轉換成二次規劃表達式再進行求解。相關數學概念參考另一篇：最優化問題基礎理論概述。

二、數學推導

1. 目標函數處理

$f(x_1, x_2)=(x_1-1)^2+(x_2-2)^2 =x_1^2+x_2^2-2x_1-4x_2+C$

其中，常數項用 $C$ 表示；

令， $X=\left[ \begin{matrix} x_1 \\ x_2 \end{matrix} \right]$ ，則
$\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}$

其中， $P=\left[\begin{matrix} 2 &0 \\ 0&2 \end{matrix} \right]，\ Q=\left[\begin{matrix} -2 \\ -4 \end{matrix} \right]$
?
關于為什么要寫成 $\frac{1}{2} X^TPX$ 形式，因為此時 $P$ 為目標函數的海塞矩陣，具體參看此鏈接。

2. 約束條件處理

$\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}$
令 $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]$ ，整理得約束條件如下：
$L_B \leqslant AX \leqslant U_B$

三、代碼編寫

??由步驟 二、數學推導 得到5個矩陣：

$P$ : 二次型矩陣(實對稱矩陣)；
$Q$ : 一次項矩陣；
$U_B$ : 上邊界矩陣；
$L_B$ : 下邊界矩陣；
$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$ 時目標函數取得最小。

本文來自互聯網用戶投稿，該文觀點僅代表作者本人，不代表本站立場。本站僅提供信息存儲空間服務，不擁有所有權，不承擔相關法律責任。
如若轉載，請注明出處：http://www.pswp.cn/pingmian/84159.shtml
繁體地址，請注明出處：http://hk.pswp.cn/pingmian/84159.shtml
英文地址，請注明出處：http://en.pswp.cn/pingmian/84159.shtml

如若內容造成侵權/違法違規/事實不符，請聯系多彩編程網進行投訴反饋email:809451989@qq.com，一經查實，立即刪除！