OSQP二次规划求解库使用说明
贺志国 2023.5.10
1. 凸二次规划的一般表达式
m
i
n
1
2
x
T
P
x
+
q
T
x
s
.
t
.
l
≤
A
x
≤
u
min \quad \frac{1}{2}x^T Px + q^Tx \qquad s.t. \quad l \leq Ax \leq u
min 2 1 x T P x + q T x s . t . l ≤ A x ≤ u 其中,
P
P
P 称为内核矩阵,要求是半正定矩阵 。正定矩阵在复数域上是共轭对称矩阵(Hermite矩阵),在实数域上是对称矩阵。因为我们在实数域内求解,
P
P
P 也是一个对称矩阵。
q
q
q 称为偏移向量,
A
A
A 称为仿射矩阵。
凸二次规划的求解算法主要有如下几种:
有效集法:Active Set Method
内点法:Interior Point Method
一阶方法:First Order Method
交替方向乘子法:Alternating Direction Method of Multipliers(ADMM)
OSQP
利用交替方向乘子法(ADMM)求解,qpOASES
则利用有效集法,从网上给出的测评报告看,OSQP
库的求解效率比qpOASES
库要高很多。
2. 使用OSQP求解凸二次规划的简单示例
为了形象地说明使用OSQP求解凸二次规划的过程,下面给出一个简单的示例:
m
i
n
1
2
x
T
[
4
1
1
2
]
x
+
[
1
1
]
T
x
min \quad \frac{1}{2}x^T\left[\begin{matrix} 4 & 1\\ 1 & 2\end{matrix}\right]x + \left[\begin{matrix} 1 \\ 1\end{matrix}\right]^Tx
min 2 1 x T [ 4 1 1 2 ] x + [ 1 1 ] T x s.t. (subject to):
[
1
0
0
]
≤
[
1
1
1
0
0
1
]
x
≤
[
1
0.7
0.7
]
\left[\begin{matrix}1 \\ 0 \\ 0 \end{matrix}\right] \leq \left[\begin{matrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{matrix}\right]x \leq \left[\begin{matrix}1 \\ 0.7 \\ 0.7 \end{matrix}\right]
1 0 0
≤
1 1 0 1 0 1
x ≤
1 0.7 0.7
本示例中,内核矩阵
P
=
[
4
1
1
2
]
P=\left[\begin{matrix} 4 & 1\\1 & 2\end{matrix}\right]
P = [ 4 1 1 2 ] 的元素全部为实数,因此它一定是对称矩阵。对于对称矩阵,只需要存储上三角矩阵即可,这也是OSQP
存储稀疏矩阵的常见做法,当然也是该库提高计算效率的有效措施。
具体而言,OSQP
库使用csc_matrix(indices, indptr, data , csc = Compressed Sparse Column)的方式来存储一个矩阵。csc_matrix存储格式定义如下:
indices is array of row indices
data is array of corresponding nonzero values
indptr points to column starts in indices and data
下面仍然以内核矩阵
P
=
[
4
1
1
2
]
P=\left[\begin{matrix} 4 & 1\\1 & 2\end{matrix}\right]
P = [ 4 1 1 2 ] 为例进行说明。首先,将内核矩阵简化为上三角矩阵:
P
=
[
4
1
0
2
]
P=\left[\begin{matrix} 4 & 1\\0 & 2\end{matrix}\right]
P = [ 4 0 1 2 ]
P
P
P 的非零元素个数为3,记为:P_nnz = 3 (nnz: number of non-zero elements)
P
P
P 中非零元素为:4, 1, 2,记为:P_x = {4, 1, 2};
P
P
P 中非零元素行索引为:0, 0, 1,这是按照先第一列,再第二列。。。,最后第N列的顺序标记。第一列的非零元素是4,行索引为0;第二列的非零元素为1, 2,行索引分别为:0, 1,于是记为:P_i = {0, 0, 1};
P
P
P 中第一列非零元素是1个,第二列非零元素是2个,记为:P_p = {0, 1, 3}。这个记法有些反人类,也很难形象理解,容我再详细解释一番。P_p元素个数为矩阵列数加1 ,P_p的第一个元素为0 ,P_p中前后两个元素的差值表示矩阵P相应列的非零元素个数 。例如:P_p = {0, 1, 3}表示:1 - 0 = 1 代表矩阵P第一列元素非零个数为1个,3 - 1 = 2 代表矩阵P第二列元素非零个数为2个。助记方法:P_i (i代表indices), P_p(p代表indptr)。
根据上述介绍,下面给出该示例的求解代码:
# include "osqp.h"
int main ( int argc, char * * argv) {
// P矩阵是对称矩阵,只存储上三角部分,下三角部分视为零值
c_float P_x[ 3 ] = { 4.0 , 1.0 , 2.0 , } ;
// P矩阵非零元素个数为3
// nnz = number of non-zero elements.
c_int P_nnz = 3 ;
// P矩阵非零元素行索引,按照先第一列,再第二列。。。
// 的顺序排列。下例表示非零元素的行序号为0、0、1
c_int P_i[ 3 ] = { 0 , 0 , 1 , } ;
// 1-0=1表示第0列非零元素个数
// 3-1=2表示第1列非零元素个数
c_int P_p[ 3 ] = { 0 , 1 , 3 , } ;
c_float q[ 2 ] = { 1.0 , 1.0 , } ;
// A矩阵非零元素
c_float A_x[ 4 ] = { 1.0 , 1.0 , 1.0 , 1.0 , } ;
// A矩阵非零元素个数为4
c_int A_nnz = 4 ;
// A矩阵非零元素的行序号为0、1、0、2(按列排序)
c_int A_i[ 4 ] = { 0 , 1 , 0 , 2 , } ;
// 2-0=2表示第0列2个元素非零
// 4-2=2表示第1列2个元素非零
c_int A_p[ 3 ] = { 0 , 2 , 4 , } ;
c_float l[ 3 ] = { 1.0 , 0.0 , 0.0 , } ;
c_float u[ 3 ] = { 1.0 , 0.7 , 0.7 , } ;
c_int n = 2 ;
c_int m = 3 ;
// Exitflag
c_int exitflag = 0 ;
// Workspace structures
OSQPWorkspace * work;
OSQPSettings * settings = ( OSQPSettings * ) c_malloc ( sizeof ( OSQPSettings) ) ;
OSQPData * data = ( OSQPData * ) c_malloc ( sizeof ( OSQPData) ) ;
// Populate data
if ( data) {
data-> n = n;
data-> m = m;
data-> P = csc_matrix ( data-> n, data-> n, P_nnz, P_x, P_i, P_p) ;
data-> q = q;
data-> A = csc_matrix ( data-> m, data-> n, A_nnz, A_x, A_i, A_p) ;
data-> l = l;
data-> u = u;
}
// Define solver settings as default
if ( settings) {
osqp_set_default_settings ( settings) ;
settings-> alpha = 1.0 ; // Change alpha parameter
}
// Setup workspace
exitflag = osqp_setup ( & work, data, settings) ;
// Solve Problem
osqp_solve ( work) ;
// Cleanup
if ( data) {
if ( data-> A) c_free ( data-> A) ;
if ( data-> P) c_free ( data-> P) ;
c_free ( data) ;
}
if ( settings) c_free ( settings) ;
return exitflag;
}