5.5.3. C++ 的MIQP建模和优化

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

首先,引入头文件:

27#include <vector>

并创建优化模型:

36    MDOEnv env = MDOEnv();
37    MDOModel model = MDOModel(env);

接下来,我们通过 MDOModel::set() 将目标函数设置为 最小化,并调用 MDOModel::addVar() 来添加四个优化变量(有关模型属性内容及其设置可参考 属性, 其他API请参考 C++ API):

44        /* Change to minimization problem. */
45        model.set(MDO_IntAttr_ModelSense, MDO_MINIMIZE);
46
47        /* Add variables. */
48        std::vector<MDOVar> x;
49        x.push_back(model.addVar(0.0, 10.0,         0.0, MDO_INTEGER, "x0"));
50        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, MDO_CONTINUOUS, "x1"));

接着,我们开始添加线性约束:

52        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, MDO_CONTINUOUS, "x3"));
53
54        /* Add constraints. */

然后,我们创建一个二次表达式 MDOQuadExpr, 再调用 MDOQuadExpr::addTerms 来设置目标函数线性部分。 obj_idx 表示线性部分的索引,obj_val 表示与 obj_idx 中的索引相对应的非零系数值,obj_nnz 代表线性部分的非零元的个数。

56        model.addConstr(1.0 * x[0] - 1.0 * x[2] + 6.0 * x[3], MDO_EQUAL, 1.0, "c1");
57        
58        /*Create a QuadExpr. */
59        MDOQuadExpr obj = MDOQuadExpr(0.0);
60
61        /* Add objective linear term.*/
62        const MDOVar obj_idx[] = { x[0], x[1], x[2], x[3]};
63        const double obj_val[] = { 1.0, 1.0, 1.0, 1.0};

然后,调用 MDOQuadExpr::addTerms 来设置目标的二次项系数 \(Q\)。 其中,qo_values 表示要添加的二次项的系数,qo_col1qo_col2 表示与qo_values相对应的二次项的第一个变量和第二个变量,qo_nnz 表示二次项中的非零元个数。

64        obj.addTerms(obj_val, obj_idx, 4);
65
66        /* Add quadratic objective matrix Q.
67         *
68         *  Note.
69         *  1. The objective function is defined as c^Tx + 1/2 x^TQx, where Q is stored with coordinate format.
70         *  2. Q will be scaled by 1/2 internally.

最后,我们调用 MDOModel::setObjective 设定优化目标与方向。

72         *

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

78

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

  1/**
  2 *  Description
  3 *  -----------
  4 *
  5 *  Linear 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 *  Subject To
 14 *   c0 : 1 x0 + 1 x1 + 2 x2 + 3 x3 >= 1
 15 *   c1 : 1 x0 - 1 x2 + 6 x3 = 1
 16 *  Bounds
 17 *    0 <= x0 <= 10
 18 *    0 <= x1
 19 *    0 <= x2
 20 *    0 <= x3
 21 *  Integers
 22 *  x0 
 23 *  End
 24 */
 25#include <iostream>
 26#include "MindoptCpp.h"
 27#include <vector>
 28
 29using namespace std;
 30
 31int main(void)
 32{
 33    /*------------------------------------------------------------------*/
 34    /* Step 1. Create environment and model.                            */
 35    /*------------------------------------------------------------------*/
 36    MDOEnv env = MDOEnv();
 37    MDOModel model = MDOModel(env);
 38    
 39    try 
 40    {
 41        /*------------------------------------------------------------------*/
 42        /* Step 2. Input model.                                             */
 43        /*------------------------------------------------------------------*/
 44        /* Change to minimization problem. */
 45        model.set(MDO_IntAttr_ModelSense, MDO_MINIMIZE);
 46
 47        /* Add variables. */
 48        std::vector<MDOVar> x;
 49        x.push_back(model.addVar(0.0, 10.0,         0.0, MDO_INTEGER, "x0"));
 50        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, MDO_CONTINUOUS, "x1"));
 51        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, MDO_CONTINUOUS, "x2"));
 52        x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, MDO_CONTINUOUS, "x3"));
 53
 54        /* Add constraints. */
 55        model.addConstr(1.0 * x[0] + 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3], MDO_GREATER_EQUAL, 1.0, "c0");
 56        model.addConstr(1.0 * x[0] - 1.0 * x[2] + 6.0 * x[3], MDO_EQUAL, 1.0, "c1");
 57        
 58        /*Create a QuadExpr. */
 59        MDOQuadExpr obj = MDOQuadExpr(0.0);
 60
 61        /* Add objective linear term.*/
 62        const MDOVar obj_idx[] = { x[0], x[1], x[2], x[3]};
 63        const double obj_val[] = { 1.0, 1.0, 1.0, 1.0};
 64        obj.addTerms(obj_val, obj_idx, 4);
 65
 66        /* Add quadratic objective matrix Q.
 67         *
 68         *  Note.
 69         *  1. The objective function is defined as c^Tx + 1/2 x^TQx, where Q is stored with coordinate format.
 70         *  2. Q will be scaled by 1/2 internally.
 71         *  3. To ensure the symmetricity of Q, user needs to input only the lower triangular part.
 72         *
 73         * Q = [ 1.0  0.5  0    0   ]
 74         *     [ 0.5  1.0  0    0   ]
 75         *     [ 0.0  0.0  1.0  0   ]
 76         *     [ 0    0    0    1.0 ]
 77         */
 78
 79        const double qo_values[] =
 80        {
 81            1.0,
 82            0.5, 1.0,
 83                    1.0, 
 84                        1.0
 85        };
 86        const MDOVar qo_col1[] = 
 87        {
 88            x[0], 
 89            x[1],   x[1],
 90                    x[2],
 91                           x[3]  
 92        };
 93        const MDOVar qo_col2[] =
 94        {
 95            x[0],
 96            x[0],   x[1],
 97                      x[2],
 98                           x[3]
 99        };
100
101        obj.addTerms(qo_values, qo_col1, qo_col2, 5);
102
103        model.setObjective(obj, MDO_MINIMIZE);
104
105        /*------------------------------------------------------------------*/
106        /* Step 3. Solve the problem.                                       */
107        /*------------------------------------------------------------------*/
108        /* Solve the problem. */
109        model.optimize();
110
111        /*------------------------------------------------------------------*/
112        /* Step 4. Retrive model status and objective.                       */
113        /* For MIP(MILP,MIQP, MIQCP) problems, if the solving process       */
114        /* terminates early due to reasons such as timeout or interruption, */
115        /* the model status will indicate termination by timeout (or        */
116        /* interruption, etc.). However, suboptimal solutions may still     */
117        /* exist, making it necessary to check the SolCount property.       */
118        /*------------------------------------------------------------------*/
119        if (model.get(MDO_IntAttr_Status) == MDO_OPTIMAL || model.get(MDO_IntAttr_Status) == MDO_SUB_OPTIMAL ||
120            model.get(MDO_IntAttr_SolCount) != 0)
121        {
122            cout << "Optimal objective value is: " << model.get(MDO_DoubleAttr_ObjVal) << endl;
123            cout << "Decision variables:" << endl;
124            int i = 0;
125            for (auto v : x)
126            {
127                cout << "x[" << i++ << "] = " << v.get(MDO_DoubleAttr_X) << endl;
128            }
129        }
130        else
131        {
132            cout<< "No feasible solution." << endl;
133        }
134        
135    }
136    catch (MDOException& e) 
137    { 
138        std::cout << "Error code = " << e.getErrorCode() << std::endl;
139        std::cout << e.getMessage() << std::endl;
140    } 
141    catch (...) 
142    { 
143        std::cout << "Error during optimization." << std::endl;
144    }
145    
146    return static_cast<int>(MDO_OKAY);
147}