5.6.8. C++ 的 MISOCP 建模和优化

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

首先,引入头文件:

26#include <iostream>

并创建优化环境和模型:

32int main(void)
33{
34    /*------------------------------------------------------------------*/
35    /* Step 1. Create environment and model.                            */
36    /*------------------------------------------------------------------*/
37    MDOEnv env = MDOEnv();

备注

MDOModel::addVar() 中,将参数 vtype 的值设置为 MDO_INTEGER 代表这个变量是整形变量。

接下来,我们通过 MDOModel::set() 将求解方向设置为 最小化。然后,调用 MDOModel::addVar() 添加五个优化变量,并指定它们的下界、上界、目标函数系数、类型和名称。

43        /* Step 2. Input model.                                             */
44        /*------------------------------------------------------------------*/
45        /* Change to minimization problem. */
46        model.set(MDO_IntAttr_ModelSense, MDO_MINIMIZE);
47
48        /* Add variables. */
49        std::vector<MDOVar> x;
50        x.push_back(model.addVar(0.0, 10.0,         1.0, MDO_INTEGER, "x0"));
51        x.push_back(model.addVar(0.0, MDO_INFINITY, 2.0, MDO_INTEGER, "x1"));
52        x.push_back(model.addVar(0.0, MDO_INFINITY, 1.0, MDO_INTEGER, "x2"));
53        x.push_back(model.addVar(0.0, MDO_INFINITY, 1.0, MDO_CONTINUOUS, "x3"));

备注

由于线性目标系数(1, 2, 1, 1, 0.5)在创建变量时已通过 addVar() 的第三个参数直接传入,因此 无需单独调用 setObjective(),且目标函数是纯线性的。

现在,我们添加线性约束。

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

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

这些约束通过重载的算术运算符来构建线性表达式并添加:

61        /* add c0 */
62        model.addConstr(1.0 * x[0] + 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3], MDO_GREATER_EQUAL, 1.0, "c0");
63
64        /* add c1 */

最后,我们添加二次(二阶锥)约束:

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

该约束使用 MDOQuadExpr 建模为二次约束。我们通过指定变量对和系数来定义二次项:

67        /* second order cone constraint */
68        int    qc2_nnz      = 3;
69        MDOVar qc2_col1[]   = { x[1], x[2], x[4]};
70        MDOVar qc2_col2[]   = { x[1], x[2], x[4]};

然后,将该二次表达式作为约束添加到模型中:

74        /* Add the second order cone constraint. c2*/
75        MDOQuadExpr c2 = MDOQuadExpr(0.0);
76        c2.addTerms(qc2_values, qc2_col1, qc2_col2, qc2_nnz);
77        model.addQConstr(c2, MDO_LESS_EQUAL, 0.0, "c2");

模型构建完成后,我们调用以下函数来求解问题:

81        /*------------------------------------------------------------------*/
82        /* Solve the problem. */

求解后,我们检查求解状态,并获取最优目标值和变量取值:

 92        /*------------------------------------------------------------------*/
 93        if (model.get(MDO_IntAttr_Status) == MDO_OPTIMAL || model.get(MDO_IntAttr_Status) == MDO_SUB_OPTIMAL ||
 94            model.get(MDO_IntAttr_SolCount) != 0)
 95        {
 96            cout << "Optimal objective value is: " << model.get(MDO_DoubleAttr_ObjVal) << endl;
 97            cout << "Decision variables:" << endl;
 98            int i = 0;
 99            for (auto v : x)
100            {
101                cout << "x[" << i++ << "] = " << v.get(MDO_DoubleAttr_X) << endl;
102            }
103        }
104        else
105        {
106            cout<< "No feasible solution." << endl;

示例 MdoMISOCPEx1.cpp 提供了完整源代码:

  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#include <iostream>
 27#include "MindoptCpp.h"
 28#include <vector>
 29
 30using namespace std;
 31
 32int main(void)
 33{
 34    /*------------------------------------------------------------------*/
 35    /* Step 1. Create environment and model.                            */
 36    /*------------------------------------------------------------------*/
 37    MDOEnv env = MDOEnv();
 38    MDOModel model = MDOModel(env);
 39
 40    try
 41    {
 42        /*------------------------------------------------------------------*/
 43        /* Step 2. Input model.                                             */
 44        /*------------------------------------------------------------------*/
 45        /* Change to minimization problem. */
 46        model.set(MDO_IntAttr_ModelSense, MDO_MINIMIZE);
 47
 48        /* Add variables. */
 49        std::vector<MDOVar> x;
 50        x.push_back(model.addVar(0.0, 10.0,         1.0, MDO_INTEGER, "x0"));
 51        x.push_back(model.addVar(0.0, MDO_INFINITY, 2.0, MDO_INTEGER, "x1"));
 52        x.push_back(model.addVar(0.0, MDO_INFINITY, 1.0, MDO_INTEGER, "x2"));
 53        x.push_back(model.addVar(0.0, MDO_INFINITY, 1.0, MDO_CONTINUOUS, "x3"));
 54        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.5, MDO_CONTINUOUS, "x4"));
 55
 56        /* Prepare model data. */
 57        /* Linear part in the objective: 1 x0 + 2 x1 + 1 x2 + 1 x3 + 0.5 x4 */
 58        int    obj_nnz   = 5;
 59        MDOVar obj_idx[] = { x[0], x[1], x[2], x[3], x[4]};
 60        double obj_val[] = { 1.0,  2.0,  1.0,  1.0, 0.5};
 61
 62        /* add c0 */
 63        model.addConstr(1.0 * x[0] + 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3], MDO_GREATER_EQUAL, 1.0, "c0");
 64
 65        /* add c1 */
 66        model.addConstr(1.0 * x[0] - 1.0 * x[2] + 6.0 * x[3], MDO_EQUAL, 1.0, "c1");
 67
 68        /* second order cone constraint */
 69        int    qc2_nnz      = 3;
 70        MDOVar qc2_col1[]   = { x[1], x[2], x[4]};
 71        MDOVar qc2_col2[]   = { x[1], x[2], x[4]};
 72        double qc2_values[] = { 1.0, 1.0, -1.0};
 73
 74        /* Add the second order cone constraint. c2*/
 75        MDOQuadExpr c2 = MDOQuadExpr(0.0);
 76        c2.addTerms(qc2_values, qc2_col1, qc2_col2, qc2_nnz);
 77        model.addQConstr(c2, MDO_LESS_EQUAL, 0.0, "c2");
 78
 79        /*------------------------------------------------------------------*/
 80        /* Step 3. Solve the problem.                                       */
 81        /*------------------------------------------------------------------*/
 82        /* Solve the problem. */
 83        model.optimize();
 84
 85        /*------------------------------------------------------------------*/
 86        /* Step 4. Retrive model status and objective.                      */
 87        /* For MIP(MILP,MIQP, MIQCP) problems, if the solving process       */
 88        /* terminates early due to reasons such as timeout or interruption, */
 89        /* the model status will indicate termination by timeout (or        */
 90        /* interruption, etc.). However, suboptimal solutions may still     */
 91        /* exist, making it necessary to check the SolCount property.       */
 92        /*------------------------------------------------------------------*/
 93        if (model.get(MDO_IntAttr_Status) == MDO_OPTIMAL || model.get(MDO_IntAttr_Status) == MDO_SUB_OPTIMAL ||
 94            model.get(MDO_IntAttr_SolCount) != 0)
 95        {
 96            cout << "Optimal objective value is: " << model.get(MDO_DoubleAttr_ObjVal) << endl;
 97            cout << "Decision variables:" << endl;
 98            int i = 0;
 99            for (auto v : x)
100            {
101                cout << "x[" << i++ << "] = " << v.get(MDO_DoubleAttr_X) << endl;
102            }
103        }
104        else
105        {
106            cout<< "No feasible solution." << endl;
107        }
108
109    }
110    catch (MDOException& e)
111    {
112        std::cout << "Error code = " << e.getErrorCode() << std::endl;
113        std::cout << e.getMessage() << std::endl;
114    }
115    catch (...)
116    {
117        std::cout << "Error during optimization." << std::endl;
118    }
119
120    return static_cast<int>(MDO_OKAY);
121}