5.6.7. C 的 MISOCP 建模和优化

在本节中，我们将使用 MindOpt C API 来建模和求解 MIQCP示例2 中的问题。

首先，引入头文件：

26#include <stdio.h>
27#include <stdlib.h>

并创建优化环境和模型：

54    /*------------------------------------------------------------------*/
55    CHECK_RESULT(MDOemptyenv(&env));
56    CHECK_RESULT(MDOstartenv(env));

接下来，我们通过 MDOsetintattr() 将求解方向设置为 最小化。然后，调用 MDOaddvar() 添加五个优化变量。它们的下界、上界、类型、名称以及在目标函数中的线性项系数都将作为参数直接传入。

备注

在函数 MDOaddvar() 中，参数 vtype 设置为 MDO_INTEGER 代表这个变量是整形变量。

59    /*------------------------------------------------------------------*/
60    /* Step 2. Input model.                                             */
61    /*------------------------------------------------------------------*/
62    /* Change to minimization problem. */
63    CHECK_RESULT(MDOsetintattr(model, MODEL_SENSE, MDO_MINIMIZE));
64
65    /* Add variables. Objective coefficients are passed directly when creating variables. */
66    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         10.0,         MDO_INTEGER,    "x0"));
67    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 2.0, 0.0,         MDO_INFINITY, MDO_INTEGER,    "x1"));
68    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         MDO_INFINITY, MDO_INTEGER,    "x2"));
69    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         MDO_INFINITY, MDO_CONTINUOUS, "x3"));
70    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 0.5, 0.0,         MDO_INFINITY, MDO_CONTINUOUS, "x4"));

备注

由于线性目标系数（1, 2, 1, 1, 0.5）在创建变量时已通过 MDOaddvar() 的第四个参数直接传入，因此 无需单独调用函数来设置或修改目标函数。

现在，我们添加线性约束。对于 C 接口，这需要为变量索引及其对应的系数准备数组。

  • 约束 c0: \(x_0 + x_1 + 2x_2 + 3x_3 \geq 1\)

  • 约束 c1: \(x_0 - x_2 + 6x_3 = 1\)

首先，准备这些约束的数据数组：

72    /* Prepare data for constraints. */
73    /* Linear constraint c0: 1 x0 + 1 x1 + 2 x2 + 3 x3 >= 1 */
74    int    c0_nnz   = 4;
75    int    c0_idx[] = {0, 1, 2, 3};
76    double c0_val[] = {1.0, 1.0, 2.0, 3.0};
77
78    /* Linear constraint c1: 1 x0 - 1 x2 + 6 x3 = 1 */
79    int    c1_nnz   = 3;
80    int    c1_idx[] = {0, 2, 3};
81    double c1_val[] = {1.0, -1.0, 6.0};

然后，使用 MDOaddconstr() 将它们添加到模型中：

89    /* Add constraints to the model. */
90    /* Add c0 */
91    CHECK_RESULT(MDOaddconstr(model, c0_nnz, c0_idx, c0_val, MDO_GREATER_EQUAL, 1.0, "c0"));
92    /* Add c1 */
93    CHECK_RESULT(MDOaddconstr(model, c1_nnz, c1_idx, c1_val, MDO_EQUAL, 1.0, "c1"));

最后，我们添加二阶锥约束：

  • c2: \(x_1^2 + x_2^2 - x_4^2 \leq 0\)

与线性约束类似，我们为二次项的变量索引和系数值准备数组：

86    int    qc2_col2[]   = {1, 2, 4};
87    double qc2_values[] = {1.0, 1.0, -1.0};

然后，使用 MDOaddqconstr() 将二次约束添加到模型中。请注意，此约束没有线性部分，因此线性分量的参数被设置为 0 和 NULL。

83    /* Second order cone constraint c2: x_1^2 + x_2^2 - x_4^2 <= 0 */
84    int    qc2_nnz      = 3;
85    int    qc2_col1[]   = {1, 2, 4};
86    int    qc2_col2[]   = {1, 2, 4};
87    double qc2_values[] = {1.0, 1.0, -1.0};

模型构建完成后，我们调用 MDOoptimize() 来求解问题：

100    /* Solve the problem. */
101    CHECK_RESULT(MDOoptimize(model));

求解后，我们检查求解状态。对于混合整数二阶锥规划问题，若求解过程因超时或中断等原因提前终止，仍可能存在次优解，因此检查 SolCount 属性很重要。然后我们获取最优目标值和变量取值。

110    /*------------------------------------------------------------------*/
111    CHECK_RESULT(MDOgetintattr(model, STATUS, &status));
112    CHECK_RESULT(MDOgetintattr(model, SOL_COUNT, &solcount));
113    if (status == MDO_OPTIMAL || status == MDO_SUB_OPTIMAL || solcount > 0)
114    {
115        CHECK_RESULT(MDOgetdblattr(model, OBJ_VAL, &obj));
116        printf("Optimal objective value is: %f\n", obj);
117        printf("Decision variables:\n");
118        for (i = 0; i < num_vars; ++i)
119        {
120            CHECK_RESULT(MDOgetdblattrelement(model, X, i, &x));
121            printf("x[%d] = %f\n", i, x);
122        }
123    }
124    else
125    {
126        printf("No feasible solution.\n");
127    }

最后，我们调用宏定义 RELEASE_MEMORY 来释放为模型和环境分配的内存：

31/* Macro to check the return code. */
32#define RELEASE_MEMORY  \
33    MDOfreemodel(model);    \

129    /*------------------------------------------------------------------*/
130    /* Step 5. Free the model.                                          */
131    /*------------------------------------------------------------------*/
132    RELEASE_MEMORY;

示例 MdoMISOCPEx1.c 提供了完整源代码：

  1/**
  2 *  Description
  3 *  -----------
  4 *
  5 *  Formulation
  6 *  -----------
  7 *
  8  Minimize
  9 *    obj: 1 x0 + 2 x1 + 1 x2 + 1 x3 + 0.5 x_4
 10 *
 11 *
 12 *  Subject To
 13 *   c0 : 1 x0 + 1 x1 + 2 x2 + 3 x3 >= 1
 14 *   c1 : 1 x0 - 1 x2 + 6 x3  = 1
 15 *   c2 : x_1^2 + x_2^2 - x_4^2 <= 0
 16 *  Bounds
 17 *    0 <= x0 <= 10
 18 *    0 <= x1
 19 *    0 <= x2
 20 *    0 <= x3
 21 *    0 <= x4
 22 *  Integer
 23 *  x0, x1, x2
 24 *  End
 25 */
 26
 27#include <stdio.h>
 28#include <stdlib.h>
 29#include "Mindopt.h"
 30
 31/* Macro to check the return code. */
 32#define RELEASE_MEMORY  \
 33    MDOfreemodel(model);    \
 34    MDOfreeenv(env);
 35#define CHECK_RESULT(code) { int res = code; if (res != 0) { fprintf(stderr, "Bad code: %d\n", res);  RELEASE_MEMORY; return (res); } }
 36#define MODEL_NAME  "MISOCPEx1"
 37#define MODEL_SENSE "ModelSense"
 38#define SOL_COUNT   "SolCount"
 39#define STATUS      "Status"
 40#define OBJ_VAL     "ObjVal"
 41#define X           "X"
 42
 43int main(void)
 44{
 45    /* Variables. */
 46    MDOenv *env;
 47    MDOmodel *model;
 48    double obj, x;
 49    int i, solcount, status;
 50    int num_vars = 5;
 51
 52    /*------------------------------------------------------------------*/
 53    /* Step 1. Create environment and model.                            */
 54    /*------------------------------------------------------------------*/
 55    CHECK_RESULT(MDOemptyenv(&env));
 56    CHECK_RESULT(MDOstartenv(env));
 57    CHECK_RESULT(MDOnewmodel(env, &model, MODEL_NAME, 0, NULL, NULL, NULL, NULL, NULL));
 58
 59    /*------------------------------------------------------------------*/
 60    /* Step 2. Input model.                                             */
 61    /*------------------------------------------------------------------*/
 62    /* Change to minimization problem. */
 63    CHECK_RESULT(MDOsetintattr(model, MODEL_SENSE, MDO_MINIMIZE));
 64
 65    /* Add variables. Objective coefficients are passed directly when creating variables. */
 66    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         10.0,         MDO_INTEGER,    "x0"));
 67    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 2.0, 0.0,         MDO_INFINITY, MDO_INTEGER,    "x1"));
 68    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         MDO_INFINITY, MDO_INTEGER,    "x2"));
 69    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 1.0, 0.0,         MDO_INFINITY, MDO_CONTINUOUS, "x3"));
 70    CHECK_RESULT(MDOaddvar(model, 0, NULL, NULL, 0.5, 0.0,         MDO_INFINITY, MDO_CONTINUOUS, "x4"));
 71
 72    /* Prepare data for constraints. */
 73    /* Linear constraint c0: 1 x0 + 1 x1 + 2 x2 + 3 x3 >= 1 */
 74    int    c0_nnz   = 4;
 75    int    c0_idx[] = {0, 1, 2, 3};
 76    double c0_val[] = {1.0, 1.0, 2.0, 3.0};
 77
 78    /* Linear constraint c1: 1 x0 - 1 x2 + 6 x3 = 1 */
 79    int    c1_nnz   = 3;
 80    int    c1_idx[] = {0, 2, 3};
 81    double c1_val[] = {1.0, -1.0, 6.0};
 82
 83    /* Second order cone constraint c2: x_1^2 + x_2^2 - x_4^2 <= 0 */
 84    int    qc2_nnz      = 3;
 85    int    qc2_col1[]   = {1, 2, 4};
 86    int    qc2_col2[]   = {1, 2, 4};
 87    double qc2_values[] = {1.0, 1.0, -1.0};
 88
 89    /* Add constraints to the model. */
 90    /* Add c0 */
 91    CHECK_RESULT(MDOaddconstr(model, c0_nnz, c0_idx, c0_val, MDO_GREATER_EQUAL, 1.0, "c0"));
 92    /* Add c1 */
 93    CHECK_RESULT(MDOaddconstr(model, c1_nnz, c1_idx, c1_val, MDO_EQUAL, 1.0, "c1"));
 94    /* Add c2: The constraint has no linear part, so the first three arguments are 0, NULL, NULL. */
 95    CHECK_RESULT(MDOaddqconstr(model, 0, NULL, NULL, qc2_nnz, qc2_col1, qc2_col2, qc2_values, MDO_LESS_EQUAL, 0.0, "c2"));
 96
 97    /*------------------------------------------------------------------*/
 98    /* Step 3. Solve the problem.                                       */
 99    /*------------------------------------------------------------------*/
100    /* Solve the problem. */
101    CHECK_RESULT(MDOoptimize(model));
102
103    /*------------------------------------------------------------------*/
104    /* Step 4. Retrive model status and objective.                      */
105    /* For MIP(MILP,MIQP, MIQCP) problems, if the solving process       */
106    /* terminates early due to reasons such as timeout or interruption, */
107    /* the model status will indicate termination by timeout (or        */
108    /* interruption, etc.). However, suboptimal solutions may still     */
109    /* exist, making it necessary to check the SolCount property.       */
110    /*------------------------------------------------------------------*/
111    CHECK_RESULT(MDOgetintattr(model, STATUS, &status));
112    CHECK_RESULT(MDOgetintattr(model, SOL_COUNT, &solcount));
113    if (status == MDO_OPTIMAL || status == MDO_SUB_OPTIMAL || solcount > 0)
114    {
115        CHECK_RESULT(MDOgetdblattr(model, OBJ_VAL, &obj));
116        printf("Optimal objective value is: %f\n", obj);
117        printf("Decision variables:\n");
118        for (i = 0; i < num_vars; ++i)
119        {
120            CHECK_RESULT(MDOgetdblattrelement(model, X, i, &x));
121            printf("x[%d] = %f\n", i, x);
122        }
123    }
124    else
125    {
126        printf("No feasible solution.\n");
127    }
128
129    /*------------------------------------------------------------------*/
130    /* Step 5. Free the model.                                          */
131    /*------------------------------------------------------------------*/
132    RELEASE_MEMORY;
133
134    return 0;
135}