5.4.2. C的SDP建模与优化

在本节中,我们将使用 MindOptC API 来建模以及求解 半定规划问题示例 中的问题。

5.4.2.1. SDP示例1

首先,需要引入头文件:

26#include "Mindopt.h"

第一步:创建优化模型:

先建立一空的MindOpt模型。

68    /* Create an empty model. */
69    CHECK_RESULT(MDOemptyenv(&env));
70    CHECK_RESULT(MDOstartenv(env));
71    CHECK_RESULT(MDOnewmodel(env, &m, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));

第二步:SDP模型输入

接着,我们将目标函数改为 最大化

77    /* Change to maximization problem. */
78    CHECK_RESULT(MDOsetintattr(m, MODEL_SENSE, MDO_MAXIMIZE));
我们利用 MDOaddpsdvar() 来创建半正定矩阵变量,并可以同时设定目标函数中其对应的矩阵系数:
  • 第一个参数为 模型

  • 第二个参数为 矩阵变量的维数 (此处为3),

  • 第三个元素代表着 矩阵系数的非零元素个数

  • 第四个和第五个参数分别代表 矩阵系数中非零元素对应的向量化索引与取值

79    /* Add matrix variable and corresponding objective coefficients. */
80    CHECK_RESULT(MDOaddpsdvar(m, dim_mat, C_size, C_nnz_indices, C_nnz_values, NULL));
接着,我们输入约束。我们利用 MDOaddpsdconstr() 创建带半正定矩阵变量的约束:
  • 第一个参数为 模型

  • 第二个参数为该约束中 线性项的个数 (此处为0),

  • 第三个参数为该约束中 矩阵项的个数

  • 第四个参数代表 所有矩阵系数的总非零元素的个数

  • 第五个和第六个参数代表该约束中 线性项对应的索引与取值

  • 第七个参数代表该约束中 矩阵变量的索引

  • 第八个参数代表该约束中 各矩阵系数的初始索引

  • 第九个和第十个参数代表 所有矩阵系数的非零元素的向量化索引与取值

  • 第十一个参数代表该 约束的意义,MDO_EQUAL (‘=’)代表等式约束,

  • 第十二个参数代表 约束的右侧值

  • 最后一个参数代表该约束中 名称

82    /* Add Constraints. */
83    CHECK_RESULT(MDOaddpsdconstr(m, 0, num_psd_var, A_size, NULL, NULL, psd_var_indices, psd_coeff_start_indices, A_nnz_indices, A_nnz_values, MDO_EQUAL, 1.0, NULL));

第三步:求解SDP模型

模型输入后,我们接着用 MDOoptimize() 来求解问题。

88    /* Solve the problem. */
89    CHECK_RESULT(MDOoptimize(m));

第四步: 取得SDP模型的解

我们利用 MDOgetintattr() 得到求解状态,即 Status 属性;利用泛用函数 MDOgetdblattr() 来取得最优的的目标函数值,即 PrimalObjVal 属性。

90    CHECK_RESULT(MDOgetintattr(m, STATUS, &status));
91    if (status == MDO_OPTIMAL)
92    {
93        CHECK_RESULT(MDOgetdblattr(m, OBJ_VAL, &obj));
94        printf("The objective value is %f\n", obj);
95    }
96    else 
97    {
98        printf("No feasible solution exists\n");
99    }

第五步: 释放模型

我们利用 MDOfreemodel()MDOfreeenv() 来释放模型。

30#define RELEASE_MEMORY  \
31    MDOfreemodel(m);    \
32    MDOfreeenv(env);
105    RELEASE_MEMORY;

示例 mdo_sdo_ex1.c 提供了完整源代码:

  1/**
  2 *  Description
  3 *  -----------
  4 *
  5 *  Semidefinite optimization (row-wise input).
  6 *
  7 *  Formulation
  8 *  -----------
  9 *
 10 *  Maximize 
 11 *    obj: tr(C X) 
 12 *  Subject To
 13 *    c0 : tr(A X) = 1
 14 *  Bounds
 15 *    X is p.s.d.
 16 *
 17 *  Matrix
 18 *    C = [ -3  0  1 ]  A = [ 3 0 1 ]
 19 *        [  0 -2  0 ]      [ 0 4 0 ]
 20 *        [  1  0 -3 ]      [ 1 0 5 ]
 21 *  End
 22 */
 23
 24#include <stdio.h>
 25#include <stdlib.h>
 26#include "Mindopt.h"
 27
 28/* Macro to check the return code. */
 29#define CHECK_RESULT(code) { int res = code; if (res != 0) { fprintf(stderr, "Bad code: %d\n", res); exit(res); } }
 30#define RELEASE_MEMORY  \
 31    MDOfreemodel(m);    \
 32    MDOfreeenv(env);
 33#define MODEL_NAME "SDP_01"
 34#define MODEL_SENSE "ModelSense"
 35#define STATUS "Status"
 36#define OBJ_VAL "PrimalObjVal"
 37
 38int main(void)
 39{
 40    /* Variables. */
 41    MDOenv *env;
 42    MDOmodel *m;
 43    double obj;
 44    int status;
 45
 46    /* Model data. */
 47    // Variable dims.
 48    int    num_mats = 1;
 49    int    dim_mat =  3 ;                                
 50
 51    // Objective coefficients. (Only the right upper part in row-major sparse format)
 52    int    C_size = 4;
 53    int    C_nnz_indices[] =  { 0,    2,   4,    8    }; 
 54    double C_nnz_values[]  =  { -3.0, 1.0, -2.0, -3.0 }; 
 55 
 56    // Constraint coefficients. (Only the right upper part in row-major sparse format)
 57    int    num_psd_var               = 1;
 58    int    psd_var_indices[]         = { 0 };            
 59    int    psd_coeff_start_indices[] = { 0 };
 60    
 61    int    A_size = 4;
 62    int    A_nnz_indices[] =  { 0,   2,    4,    8    }; 
 63    double A_nnz_values[]  =  { 3.0, 1.0,  4.0,  5.0  }; 
 64    
 65    /*------------------------------------------------------------------*/
 66    /* Step 1. Create a model and change the parameters.                */
 67    /*------------------------------------------------------------------*/
 68    /* Create an empty model. */
 69    CHECK_RESULT(MDOemptyenv(&env));
 70    CHECK_RESULT(MDOstartenv(env));
 71    CHECK_RESULT(MDOnewmodel(env, &m, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));
 72
 73    /*------------------------------------------------------------------*/
 74    /* Step 2. Input model.                                             */
 75    /*------------------------------------------------------------------*/
 76    /* Change to maximization problem. */
 77    CHECK_RESULT(MDOsetintattr(m, MODEL_SENSE, MDO_MAXIMIZE));
 78
 79    /* Add matrix variable and corresponding objective coefficients. */
 80    CHECK_RESULT(MDOaddpsdvar(m, dim_mat, C_size, C_nnz_indices, C_nnz_values, NULL));
 81
 82    /* Add Constraints. */
 83    CHECK_RESULT(MDOaddpsdconstr(m, 0, num_psd_var, A_size, NULL, NULL, psd_var_indices, psd_coeff_start_indices, A_nnz_indices, A_nnz_values, MDO_EQUAL, 1.0, NULL));
 84  
 85    /*------------------------------------------------------------------*/
 86    /* Step 3. Solve the problem and populate the result.               */
 87    /*------------------------------------------------------------------*/
 88    /* Solve the problem. */
 89    CHECK_RESULT(MDOoptimize(m));
 90    CHECK_RESULT(MDOgetintattr(m, STATUS, &status));
 91    if (status == MDO_OPTIMAL)
 92    {
 93        CHECK_RESULT(MDOgetdblattr(m, OBJ_VAL, &obj));
 94        printf("The objective value is %f\n", obj);
 95    }
 96    else 
 97    {
 98        printf("No feasible solution exists\n");
 99    }
100
101    /*------------------------------------------------------------------*/
102    /* Step 4. Free the model.                                          */
103    /*------------------------------------------------------------------*/
104    /* Free the model. */
105    RELEASE_MEMORY;
106
107    return 0;
108}

5.4.2.2. SDP示例2

首先,需要引入头文件:

32#include "Mindopt.h"

第一步:创建MindOpt模型

我们先建立一空的MindOpt模型。

94    /* Create an empty model. */
95    CHECK_RESULT(MDOemptyenv(&env));
96    CHECK_RESULT(MDOstartenv(env));
97    CHECK_RESULT(MDOnewmodel(env, &m, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));

第二步:SDP模型输入

接着,我们将目标函数改为 最大化

102    /* Change to maximization problem. */
103    CHECK_RESULT(MDOsetintattr(m, MODEL_SENSE, MDO_MAXIMIZE));

我们利用 MDOaddvar() 创建两个线性变量,并设定他们的上下界与目标函数中的系数。

105    /* Add linear variable and corresponding objective coefficients. */
106    CHECK_RESULT(MDOaddvar(m, 0, NULL, NULL, 0.0, 0.0, MDO_INFINITY, MDO_CONTINUOUS, NULL));
107    CHECK_RESULT(MDOaddvar(m, 0, NULL, NULL, 0.0, 0.0, MDO_INFINITY, MDO_CONTINUOUS, NULL));
我们再利用 MDOaddpsdvar() 来创建两个半正定矩阵变量,并同时设定目标函数中其对应的矩阵系数:
  • 第一个参数为 模型

  • 第二个参数为 矩阵变量的维数 (此处为3),

  • 第三个元素代表着 矩阵系数的非零元素个数

  • 第四个和第五个参数分别代表 矩阵系数中非零元素对应的向量化索引与取值

109    /* Add psd variable and corresponding objective coefficients. */
110    CHECK_RESULT(MDOaddpsdvar(m, dim_mats[0], C0_size, C0_nnz_indices, C0_nnz_values, NULL));
111    CHECK_RESULT(MDOaddpsdvar(m, dim_mats[1], C1_size, C1_nnz_indices, C1_nnz_values, NULL));
接着,我们输入约束。我们利用 MDOaddpsdconstr() 创建带半正定矩阵变量的约束:
  • 第一个参数为 模型

  • 第二个参数为该约束中 线性项的个数 (此处为0),

  • 第三个参数为该约束中 矩阵项的个数

  • 第四个参数代表 所有矩阵系数的总非零元素的个数

  • 第五个和第六个参数代表该约束中 线性项对应的索引与取值

  • 第七个参数代表该约束中 矩阵变量的索引

  • 第八个参数代表该约束中 各矩阵系数的初始索引

  • 第九个和第十个参数代表 所有矩阵系数的非零元素的向量化索引与取值

  • 第十一个参数代表该 约束的意义,MDO_EQUAL(‘=’)表示等式约束,

  • 第十二个参数代表 约束的右侧值

  • 最后一个参数代表该约束中 名称

113    /* Add constraints coefficients. */
114    CHECK_RESULT(MDOaddpsdconstr(m, row0_size, num_psd_var[0], A0_size, row0_idx, row0_val, A0_psd_var_indices, A0_psd_coeff_start_indices, A0_nnz_indices, A0_nnz_values, MDO_EQUAL, 1.0, NULL));
115    CHECK_RESULT(MDOaddpsdconstr(m, row1_size, num_psd_var[1], A1_size, row1_idx, row1_val, A1_psd_var_indices, A1_psd_coeff_start_indices, A1_nnz_indices, A1_nnz_values, MDO_EQUAL, 2.0, NULL));

第三步:求解SDP模型

模型输入后,我们接着用 MDOoptimize() 来求解问题。

120    /* Solve the problem. */
121    CHECK_RESULT(MDOoptimize(m));

第四步: 取得SDP模型的解

我们利用 MDOgetintattr() 得到求解状态,即 Status 属性;利用泛用函数 MDOgetdblattr() 来取得最优的的目标函数值,即 PrimalObjVal 属性。

122    CHECK_RESULT(MDOgetintattr(m, STATUS, &status));
123    if (status == MDO_OPTIMAL)
124    {
125        CHECK_RESULT(MDOgetdblattr(m, OBJ_VAL, &obj));
126        printf("The objective value is %f\n", obj);
127    }
128    else 
129    {
130        printf("No feasible solution exists\n");
131    }

第五步: 释放模型

我们利用 MDOfreemodel()MDOfreeenv() 来释放模型。

36#define RELEASE_MEMORY  \
37    MDOfreemodel(m);    \
38    MDOfreeenv(env);
137    RELEASE_MEMORY;

示例 mdo_sdo_ex2.c 提供了完整源代码:

  1/**
  2 *  Description
  3 *  -----------
  4 *
  5 *  Semidefinite optimization (row-wise input).
  6 *
  7 *  Formulation
  8 *  -----------
  9 * 
 10 *  Maximize
 11 *    obj: tr(C0 X0)   + tr(C1 X1)    + 0 x0 + 0 x1
 12 *  Subject To
 13 *    c0 : tr(A00 X0)                +  1 x0        = 1
 14 *    c1 :               tr(A11 X1)          + 1 x1 = 2
 15 *  Bounds
 16 *    0 <= x0
 17 *    0 <= x1
 18 *    X0, X1 are p.s.d.
 19 *
 20 *  Matrix
 21 *    C0 =  [ 2 1 ]   A00 = [ 3 1 ]
 22 *          [ 1 2 ]         [ 1 3 ]
 23 *
 24 *    C1 = [ 3 0 1 ]  A11 = [ 3 0 1 ]
 25 *         [ 0 2 0 ]        [ 0 4 0 ]
 26 *         [ 1 0 3 ]        [ 1 0 5 ]
 27 *  End
 28 */
 29
 30#include <stdio.h>
 31#include <stdlib.h>
 32#include "Mindopt.h"
 33
 34/* Macro to check the return code. */
 35#define CHECK_RESULT(code) { int res = code; if (res != 0) { fprintf(stderr, "Bad code: %d\n", res); exit(res); } }
 36#define RELEASE_MEMORY  \
 37    MDOfreemodel(m);    \
 38    MDOfreeenv(env);
 39#define MODEL_NAME "SDP_02"
 40#define MODEL_SENSE "ModelSense"
 41#define STATUS "Status"
 42#define OBJ_VAL "PrimalObjVal"
 43
 44int main(void)
 45{
 46    /* Variables. */
 47    MDOenv *env;
 48    MDOmodel *m;
 49    double obj;
 50    int status;
 51
 52    /* Model data. */
 53    // linear variables.
 54    int     num_cols = 2;    
 55
 56    // psd variables.
 57    int    num_mats    = 2;
 58    int    dim_mats[]  = { 2, 3 };
 59    
 60    // Objective coefficients.  (Only the right upper part in row-major sparse format)
 61    int    C0_size = 3;
 62    int    C0_nnz_indices[] =  { 0,    1,   3   }; 
 63    double C0_nnz_values[]  =  { 2.0,  1.0, 2.0 }; 
 64
 65    int    C1_size = 4;
 66    int    C1_nnz_indices[] =  { 0,    2,   4,    8   }; 
 67    double C1_nnz_values[]  =  { 3.0,  1.0, 2.0,  3.0 };  
 68    
 69    // Constraints (linear).
 70    int    row0_size  = 1;
 71    int    row0_idx[] = { 0 };
 72    double row0_val[] = { 1.0 };
 73    int    row1_size  = 1;
 74    int    row1_idx[] = { 1 };
 75    double row1_val[] = { 1.0 };
 76
 77    // Constraints (psd).  (Only the right upper part in row-major sparse format)
 78    int    num_psd_var[]                = { 1, 1 };
 79    int    A0_psd_var_indices[]         = { 0 }; 
 80    int    A0_psd_coeff_start_indices[] = { 0 };
 81    int    A0_size                      = 3;
 82    int    A0_nnz_indices[]             = { 0,   1,    3   }; 
 83    double A0_nnz_values[]              = { 3.0, 1.0,  3.0 };
 84
 85    int    A1_psd_var_indices[]         = { 1 };
 86    int    A1_psd_coeff_start_indices[] = { 0 }; 
 87    int    A1_size                      = 4;
 88    int    A1_nnz_indices[]             = { 0,   2,    4,    8   };
 89    double A1_nnz_values[]              = { 3.0, 1.0,  4.0,  5.0 }; 
 90 
 91    /*------------------------------------------------------------------*/
 92    /* Step 1. Create a model and change the parameters.                */
 93    /*------------------------------------------------------------------*/
 94    /* Create an empty model. */
 95    CHECK_RESULT(MDOemptyenv(&env));
 96    CHECK_RESULT(MDOstartenv(env));
 97    CHECK_RESULT(MDOnewmodel(env, &m, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));
 98    
 99    /*------------------------------------------------------------------*/
100    /* Step 2. Input model.                                             */
101    /*------------------------------------------------------------------*/
102    /* Change to maximization problem. */
103    CHECK_RESULT(MDOsetintattr(m, MODEL_SENSE, MDO_MAXIMIZE));
104    
105    /* Add linear variable and corresponding objective coefficients. */
106    CHECK_RESULT(MDOaddvar(m, 0, NULL, NULL, 0.0, 0.0, MDO_INFINITY, MDO_CONTINUOUS, NULL));
107    CHECK_RESULT(MDOaddvar(m, 0, NULL, NULL, 0.0, 0.0, MDO_INFINITY, MDO_CONTINUOUS, NULL));
108
109    /* Add psd variable and corresponding objective coefficients. */
110    CHECK_RESULT(MDOaddpsdvar(m, dim_mats[0], C0_size, C0_nnz_indices, C0_nnz_values, NULL));
111    CHECK_RESULT(MDOaddpsdvar(m, dim_mats[1], C1_size, C1_nnz_indices, C1_nnz_values, NULL));
112
113    /* Add constraints coefficients. */
114    CHECK_RESULT(MDOaddpsdconstr(m, row0_size, num_psd_var[0], A0_size, row0_idx, row0_val, A0_psd_var_indices, A0_psd_coeff_start_indices, A0_nnz_indices, A0_nnz_values, MDO_EQUAL, 1.0, NULL));
115    CHECK_RESULT(MDOaddpsdconstr(m, row1_size, num_psd_var[1], A1_size, row1_idx, row1_val, A1_psd_var_indices, A1_psd_coeff_start_indices, A1_nnz_indices, A1_nnz_values, MDO_EQUAL, 2.0, NULL));
116  
117    /*------------------------------------------------------------------*/
118    /* Step 3. Solve the problem and populate the result.               */
119    /*------------------------------------------------------------------*/
120    /* Solve the problem. */
121    CHECK_RESULT(MDOoptimize(m));
122    CHECK_RESULT(MDOgetintattr(m, STATUS, &status));
123    if (status == MDO_OPTIMAL)
124    {
125        CHECK_RESULT(MDOgetdblattr(m, OBJ_VAL, &obj));
126        printf("The objective value is %f\n", obj);
127    }
128    else 
129    {
130        printf("No feasible solution exists\n");
131    }
132
133    /*------------------------------------------------------------------*/
134    /* Step 4. Free the model.                                          */
135    /*------------------------------------------------------------------*/
136    /* Free the model. */
137    RELEASE_MEMORY;
138
139    return 0;
140}