5.4.2. C 的QCP建模和优化

在本节中,我们将使用 MindOpt C API,以按行输入的形式来建模以及求解 二次约束规划问题示例 中的问题。

首先,引入头文件:

27#include "Mindopt.h"

并创建优化模型:

78    CHECK_RESULT(MDOemptyenv(&env));
79    CHECK_RESULT(MDOstartenv(env));
80    CHECK_RESULT(MDOnewmodel(env, &model, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));

接下来,我们通过 MDOsetIntAttr() 将目标函数设置为 最小化,并调用 MDOaddvar() 来添加四个优化变量,定义其下界、上界、名称和类型,以及其在目标函数中线性项的系数(关于 MDOsetIntAttr()MDOaddvar() 的详细使用方式,请参考 属性):

86    /* Change to minimization problem. */
87    CHECK_RESULT(MDOsetintattr(model, MODEL_SENSE, MDO_MINIMIZE));
88
89    /* Add variables. */
90    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0,         10.0, MDO_CONTINUOUS, "x0"));
91    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x1"));
92    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x2"));
93    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x3"));

Note

矩阵的非零元随后将按 输入;因此,MDOaddvar() 中,与约束矩阵相关联的参数 sizeindicesvalue 分别用 0NULLNULL 代替(换句话说,此时问题无约束)。

接下来我们将添加二次规划中的目标函数的二次项系数。我们使用以下三列数组来定义:其中 qo_col1qo_col2 分别记录二次项中所有非零项的两个变量索引,而 qo_values 是与之相对应的非零系数值。

49    /* Quadratic part in objective: 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1] */
50    int    qo_nnz      = 5;
51    int    qo_col1[]   = { 0,   1,   2,   3,   0 };
52    int    qo_col2[]   = { 0,   1,   2,   3,   1 };
53    double qo_values[] = { 0.5, 0.5, 0.5, 0.5, 0.5 };

我们调用 MDOaddqpterms() 设置目标的二次项:

95    /* Add quadratic objective term. */
96    CHECK_RESULT(MDOaddqpterms(model, qo_nnz, qo_col1, qo_col2, qo_values));

以下我们将开始添加模型中的两个二次约束。其中的线性项与添加线性约束的方式类似,我们使用以下数组来定义线性项;其中, row1_idxrow2_idx 分别表示第一和第二个约束中非零元素的位置(索引),而 row1_valrow2_val 则是与之相对应的非零数值。

55    /* Linear part in the first constraint: 1 x0 + 1 x1 + 2 x2 + 3 x3 */
56    int    row1_nnz   = 4;
57    int    row1_idx[] = { 0,   1,   2,   3   };
58    double row1_val[] = { 1.0, 1.0, 2.0, 3.0 };
65    /* Linear part in the second constraint: 1 x0 - 1 x2 + 6 x3 */
66    int    row2_nnz   = 3;
67    int    row2_idx[] = { 0,    2,   3   };
68    double row2_val[] = { 1.0, -1.0, 6.0 };

二次约束中的二次项与目标函数中的二次项类似:第一个二次约束中我们用 qc1_col1qc1_col2 分别记录二次项中所有非零项的两个变量索引,而 qc1_values 是与之相对应的非零系数值。第二个二次约束类似。

59    /* Quadratic part in the first constraint: - 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1] */
60    int    qc1_nnz      = 5;
61    int    qc1_col1[]   = { 0,    1,    2,    3,    0 };
62    int    qc1_col2[]   = { 0,    1,    2,    3,    1 };
63    double qc1_values[] = { -0.5, -0.5, -0.5, -0.5, -0.5 };
69    /* Quadratic part in the second constraint: 1/2 [x1^2] */
70    int    qc2_nnz      = 1;
71    int    qc2_col1[]   = { 1 };
72    int    qc2_col2[]   = { 1 };
73    double qc2_values[] = { 0.5 };

我们调用 MDOaddqconstr() 来添加二次约束:

 98    /* Add quadratic constraints. */
 99    CHECK_RESULT(MDOaddqconstr(model, row1_nnz, row1_idx, row1_val, qc1_nnz, qc1_col1, qc1_col2, qc1_values, MDO_GREATER_EQUAL, 1.0, "c0"));
100    CHECK_RESULT(MDOaddqconstr(model, row2_nnz, row2_idx, row2_val, qc2_nnz, qc2_col1, qc2_col2, qc2_values, MDO_LESS_EQUAL,    1.0, "c1"));

问题输入完成后,再调用 MDOoptimize() 求解优化问题:

105    /* Solve the problem. */
106    CHECK_RESULT(MDOoptimize(model));

然后,我们可以通过获取属性值的方式来获取对应的优化目标值objective和变量的取值:

108    CHECK_RESULT(MDOgetintattr(model, STATUS, &status));
109    if (status == MDO_OPTIMAL) 
110    {
111        CHECK_RESULT(MDOgetdblattr(model, OBJ_VAL, &obj));
112        printf("The optimal objective value is: %f\n", obj);
113        for (int i = 0; i < 4; ++i) 
114        {
115            CHECK_RESULT(MDOgetdblattrelement(model, X, i, &x));
116            printf("x[%d] = %f\n", i, x);
117        }
118    } 
119    else 
120    {
121        printf("No feasible solution.\n");
122    }

最后,调用 MDOfreemodel()MDOfreeenv() 来释放模型:

30#define RELEASE_MEMORY  \
31    MDOfreemodel(model);    \
32    MDOfreeenv(env);
127    RELEASE_MEMORY;

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

  1/**
  2 *  Description
  3 *  -----------
  4 *
  5 *  Quadratically constrained quadratic optimization (row-wise input).
  6 *
  7 *  Formulation
  8 *  -----------
  9 *
 10 *  Minimize
 11 *    obj: 1 x0 + 1 x1 + 1 x2 + 1 x3
 12 *         + 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1]
 13 *
 14 *  Subject To
 15 *   c0 : 1 x0 + 1 x1 + 2 x2 + 3 x3 - 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1] >= 1
 16 *   c1 : 1 x0 - 1 x2 + 6 x3 + 1/2 [x1^2] <= 1
 17 *  Bounds
 18 *    0 <= x0 <= 10
 19 *    0 <= x1
 20 *    0 <= x2
 21 *    0 <= x3
 22 *  End
 23 */
 24
 25#include <stdio.h>
 26#include <stdlib.h>
 27#include "Mindopt.h"
 28
 29/* Macro to check the return code */
 30#define RELEASE_MEMORY  \
 31    MDOfreemodel(model);    \
 32    MDOfreeenv(env);
 33#define CHECK_RESULT(code) { int res = code; if (res != 0) { fprintf(stderr, "Bad code: %d\n", res);  RELEASE_MEMORY; return (res); } }
 34#define MODEL_NAME  "QCP_01"
 35#define MODEL_SENSE "ModelSense"
 36#define STATUS      "Status"
 37#define OBJ_VAL     "ObjVal"
 38#define X           "X"
 39
 40int main(void)
 41{
 42    /* Variables. */
 43    MDOenv *env;
 44    MDOmodel *model;
 45    double obj, x;
 46    int status, i;
 47
 48    /* Prepare model data. */
 49    /* Quadratic part in objective: 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1] */
 50    int    qo_nnz      = 5;
 51    int    qo_col1[]   = { 0,   1,   2,   3,   0 };
 52    int    qo_col2[]   = { 0,   1,   2,   3,   1 };
 53    double qo_values[] = { 0.5, 0.5, 0.5, 0.5, 0.5 };
 54
 55    /* Linear part in the first constraint: 1 x0 + 1 x1 + 2 x2 + 3 x3 */
 56    int    row1_nnz   = 4;
 57    int    row1_idx[] = { 0,   1,   2,   3   };
 58    double row1_val[] = { 1.0, 1.0, 2.0, 3.0 };
 59    /* Quadratic part in the first constraint: - 1/2 [ x0^2 + x1^2 + x2^2 + x3^2 + x0 x1] */
 60    int    qc1_nnz      = 5;
 61    int    qc1_col1[]   = { 0,    1,    2,    3,    0 };
 62    int    qc1_col2[]   = { 0,    1,    2,    3,    1 };
 63    double qc1_values[] = { -0.5, -0.5, -0.5, -0.5, -0.5 };
 64
 65    /* Linear part in the second constraint: 1 x0 - 1 x2 + 6 x3 */
 66    int    row2_nnz   = 3;
 67    int    row2_idx[] = { 0,    2,   3   };
 68    double row2_val[] = { 1.0, -1.0, 6.0 };
 69    /* Quadratic part in the second constraint: 1/2 [x1^2] */
 70    int    qc2_nnz      = 1;
 71    int    qc2_col1[]   = { 1 };
 72    int    qc2_col2[]   = { 1 };
 73    double qc2_values[] = { 0.5 };
 74
 75     /*------------------------------------------------------------------*/
 76    /* Step 1. Create environment and model.                            */
 77    /*------------------------------------------------------------------*/
 78    CHECK_RESULT(MDOemptyenv(&env));
 79    CHECK_RESULT(MDOstartenv(env));
 80    CHECK_RESULT(MDOnewmodel(env, &model, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));
 81
 82
 83    /*------------------------------------------------------------------*/
 84    /* Step 2. Input model.                                             */
 85    /*------------------------------------------------------------------*/
 86    /* Change to minimization problem. */
 87    CHECK_RESULT(MDOsetintattr(model, MODEL_SENSE, MDO_MINIMIZE));
 88
 89    /* Add variables. */
 90    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0,         10.0, MDO_CONTINUOUS, "x0"));
 91    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x1"));
 92    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x2"));
 93    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0, MDO_INFINITY, MDO_CONTINUOUS, "x3"));
 94
 95    /* Add quadratic objective term. */
 96    CHECK_RESULT(MDOaddqpterms(model, qo_nnz, qo_col1, qo_col2, qo_values));
 97
 98    /* Add quadratic constraints. */
 99    CHECK_RESULT(MDOaddqconstr(model, row1_nnz, row1_idx, row1_val, qc1_nnz, qc1_col1, qc1_col2, qc1_values, MDO_GREATER_EQUAL, 1.0, "c0"));
100    CHECK_RESULT(MDOaddqconstr(model, row2_nnz, row2_idx, row2_val, qc2_nnz, qc2_col1, qc2_col2, qc2_values, MDO_LESS_EQUAL,    1.0, "c1"));
101    
102    /*------------------------------------------------------------------*/
103    /* Step 3. Solve the problem and populate optimization result.                */
104    /*------------------------------------------------------------------*/
105    /* Solve the problem. */
106    CHECK_RESULT(MDOoptimize(model));
107        
108    CHECK_RESULT(MDOgetintattr(model, STATUS, &status));
109    if (status == MDO_OPTIMAL) 
110    {
111        CHECK_RESULT(MDOgetdblattr(model, OBJ_VAL, &obj));
112        printf("The optimal objective value is: %f\n", obj);
113        for (int i = 0; i < 4; ++i) 
114        {
115            CHECK_RESULT(MDOgetdblattrelement(model, X, i, &x));
116            printf("x[%d] = %f\n", i, x);
117        }
118    } 
119    else 
120    {
121        printf("No feasible solution.\n");
122    }
123 
124    /*------------------------------------------------------------------*/
125    /* Step 4. Free the model.                                          */
126    /*------------------------------------------------------------------*/
127    RELEASE_MEMORY;
128       
129    return 0;
130}