C++嵌套调用 用Python 脚本写的 基于Gurobi 的解数学模型的经验记录

The Experience of How to Solve Mathematical Model written in Python, based on Gurobi in C++ simulation program.

一、前言

标题写得有点拗口，总共涉及到3个工具：Python, Gurobi, C++，其中 Gurobi 是一个类似于 IBM 的 CPlex 的解数学模型的工具，下文有详述。

这篇记录是把我近一段时间 用 C++ 调用 由python 脚本 写的基于 Gurobi optimizer 的 程序的经验分享出来。

我当初开始试着解决问题时，搜索了整个网络，国语写成的资料，最多涉及到2项，如 “如何用C++调用 python 脚本” 或者 “如何 调用 Gurobi 的 C++ 接口”，又或者 “ 如何调用 Gurobi 的 python 接口”。 没有一个符合我的要求的。木有办法，只好去参考原版的
Gurobi Manual 与 python 27 的 tutorial， 然后一步步地做实验，探索用法。然后结合着 Google Group （https://groups.google.com/forum/#!forum/gurobi）
的帖子，去借鉴别人的经验。虽然这个过程走了很多的弯路，最后终于把自己的问题解决了！

首先，把我的需求说明清楚。

我遇到的需求：

1）已经有一个数学model，但是不知道是否正确，所以需要用一个简单的小例子快速地验证；

2）如果 model 正确了，我需要把它封装一下，嵌到我的C++工程里，这样可以在模拟程序里一轮一轮地反复调用，记录多次调用
model 解出来的optimal 值，同时可以让模拟程序里设计的算法执行在由同样的输入数据的情况下，最终可以比较算法的结果与最优值之间的 performance 差别。

前边 之所以说，走了弯路，而且走了两次：一是因为最开始时，我没有挑选对优化工具，用了国内常用的Lingo，因为在国内网站上搜到很多诸如“谢金星Lingo优化工具讲解” ，擦，这个工具的脚本语言真是不好用，不是贬损人家，是我用着真不顺手，花了好几天的时间还用不好。渐渐生出换解model
工具的念头。请教师兄，他推荐我使用世界上优化效率最好的 Gurobi，比著名的CPlex还要好用。最开始不信，google一下再说，结果第一条便是“Gurobi —
The overall fastest and best supported solver available”，而且找到了几种主流的优化工具效率的比较的柱状图，不是重点，此处略去。

开读 Gurobi 的 manual，我发现，它的使用很灵活，提供了若干种主流语言的 interfaces， 如 C, C++, Java, Python, .net等等，比较适合嵌套在用户的程序里去被调用，来解决优化问题。

紧接着，我又走了一次弯路，因为我最熟悉C++，然后就参照着 Gurobi 提供的C++ examples，与 quick start Refs，用了C++ interfaces 来写解 model 的程序。经过千辛万苦，终于把model的solving
程序封装并写好了调用。一试用效果，小规模的程序很快，model 也被验证了正确性，心里窃喜。可是一旦上了规模，结果程序慢得像蜗牛，因为输入的数据上了规模后，用来被测试的一个Model 里的变量与 constraints 的数量呈非线性增长，有时 constraints 的数量甚至达到了数百万个，解起来很费时，根本没法用在模拟程序里进行成百上千次的反复调用。究其原因，可能性是1）我的水平不够，写的C++程序效率奇差；2）我的数学model本身的特殊性造成的，当变量规模大了，constraints就异常多；3）也可能是因为
Gurobi 所提供的C++ interfaces 在解我的 model 时，效果不好。

木有办法，再次请教师兄，他建议我使用 python试试（现在看来，真是前辈的经验宝贵啊）。于是，就现学现用python了。对照着师兄给的sample，慢慢将我的model写成用python脚本来解。又经过1天的埋头coding，终于写好了单独的一个用
python 脚本写好的 model solving 程序。开试效果，哇，只要规模不是太大，几乎瞬间被解出。怎么回事？ 现在看来刚才用 C++ 调用效果不好的3个可能的原因，应该是第3条被验证了。

好，那我就用 python 脚本写成的 model-solving 程序来解最优值。

Then, the next problem occurs: How to use my C++ project to call this python model-solving script ?

到了这一步，这个问题就简单多了，只需要解决用C++来调用 python 脚本就 over了。很快，这个问题被我解决了。

以上便是我的此次的问题解决路线，弯路有时候必须走一走，到终点时，才知道它原来是一条弯路。尽管觉得费了很大力气，做了很多无用功，但是，我倒不觉得这是在浪费时间与精力。弯路才能让人领略到更多风景。

Ok, 废话不多说了，我把我的 用 python 写的 gurobi 脚本 与 一部分 C++ 调用这个脚本的关键程序贴在最后，如果对谁有帮助，那我的汗水也会增加一分重量。

贴出来的程序分为2个主要部分：1）python 写的 model solving 脚本示例； 2）用C++调用 Python 脚本的一些 api 的用法与封装示例.

关于用到的数学模型，如下图的 optimization 问题：



其中，各个限制条件请具体参考论文“Optimal VM Placement in Data Centers with Architectural and Resource Constraints”。这个优化问题有点类似于背包问题或者Bin
packing 问题，又或者minimum k-cut problem。解它的脚本原样贴出。

2015-Otc-03 补充: 另外一个更好的解模的例子的 Source-code 与相应的论文，请参考如下publication (Title: The_Source_Code_sample_to_solve_an_optimization_model_using_Gurobi )
URL: https://www.researchgate.net/publication/282362692
// ================================================================
//

二、程序示例与解释

1) The implementation of mathematical model based on Gurobi, in Python
script:

简单解释：此Python脚本主要提供了一个函数 SolveModel
( Input_Arg )，把 解model 所需的输入数据编成字符串当做 SolveModel function 的参数输入， 如 #Input_Arg = '3*2*3*1*0,2,4;2,0,3;3,4,0*5,5*2,2,3*1,0,0;1,0,0'，然后，函数内部会解析这个字符串，获得解模所需的初始化数据。解完后，如果需要返回所有的具体结果，那么可以将
optimal 的结果存入一个 dict 结构里，传出，然后用相关的 c++ 嵌套 python 的专用 api 来解析返回的 dict，从而在C++工程里获取最优解。这里，我只是返回一个 optimal value 就行了，所以解析返回结果的时候也就简单一些，否则如果返回一个 dict，那就需要按照你自己定义的格式来解析 dict 里的内容了。

from gurobipy import *
Cost_sameS = 1          # Cs in paper
Cost_diffS_sameK = 3    # Ck in paper
Cost_diffK = 5          # Ci in paper

def matrixMul(A, B):  # A function that gets the multiplication result between two matrices.
return [[sum(a * b for a, b in zip(a, b)) for b in zip(*B)] for a in A]

#Input_Arg = '3*2*3*1*0,2,4;2,0,3;3,4,0*5,5*2,2,3*1,0,0;1,0,0'

def SolveModel( Input_Arg ):
##    print '\n------- Input_Arg:'
##    print Input_Arg
##    print '----------------'
# ============================================ I. Load the argument ============================== :
Dict_Input_strArgs = {}
Dict_Input_strArgs = Input_Arg.strip().split('*')

CNT_K = int( Dict_Input_strArgs[0] );         # -------------- 1 load the CNT_K
CNT_S = int( Dict_Input_strArgs[1] );         # -------------- 2 load the CNT_S
CNT_V = int( Dict_Input_strArgs[2] );         # -------------- 3 load the CNT_V
CNT_R = int( Dict_Input_strArgs[3] );         # -------------- 4 load the CNT_R

# -------------- 5 load the  matT_ij
if Dict_Input_strArgs[4] is not False:
matT_ij = Dict_Input_strArgs[4].strip().split(';')
# -------------- 6 load the  matCap_rs
if Dict_Input_strArgs[5] is not False:
matCap_rs = Dict_Input_strArgs[5].strip().split(';')
# -------------- 7 load the  matQ_rv
if Dict_Input_strArgs[6] is not False:
matQ_rv = Dict_Input_strArgs[6].strip().split(';')
# -------------- 8 load the  matM_sk
if Dict_Input_strArgs[7] is not False:
matM_sk = Dict_Input_strArgs[7].strip().split(';')

# =========================================== II. Analyze the argument =============================== :
# --- T_ij
T_ij = []  # define the matrix__T_ij
for item in matT_ij:
Tij_row = item.strip().split(',')
intTij_row = []
for i in range(len(Tij_row)):
intTij_row.append( int( Tij_row[i] ) )
T_ij.append( intTij_row )

# --- Cap_rs
Cap_rs = []
for item in matCap_rs:
Cap_rs_row = item.strip().split(',')
intCap_rs_row = []
for i in range(len(Cap_rs_row)):
intCap_rs_row.append( int( Cap_rs_row[i]) )
Cap_rs.append( intCap_rs_row )

# --- Q_rv
Q_rv = []
for item in matQ_rv:
Q_rv_row = item.strip().split(',')
intQ_rv_row = []
for i in range(len(Q_rv_row)):
intQ_rv_row.append( int( Q_rv_row[i]) )
Q_rv.append( intQ_rv_row )

# --- M_sk
M_sk = []
for item in matM_sk:
M_sk_row = item.strip().split(',')
intM_sk_row = []
for i in range(len(M_sk_row)):
intM_sk_row.append( int( M_sk_row[i]) )
M_sk.append( intM_sk_row )

##    # Given values
##    T_ij = [ [0,2,3], [2,0,4], [3,4,0] ]		# Traffic_amount between VM_i and VM_j
##    Cap_rs = [[5,5]]    # Cap_rs is the capacity of Resource_r on Server_s
##    Q_rv = [[2,2,3]]  # Resource_r required by VM_v
##
##    # ---- Constant-matrixs
##    M_sk = [[1,0,0], [1,0,0]]
##    print '------------------M_sk:'
##    print M_sk
##    print '-------------------\n'

# ==========================================================================================

# ======================== 1 Create optimization model.
m = Model('VMs_placement')
m.setParam( 'OutputFlag', False )

# ======================== 2 Add vars into model.
dM_vs = {}
for v in range(CNT_V):
for s in range(CNT_S):
dM_vs[v,s] = m.addVar( vtype=GRB.BINARY, name='dM_vs_%s_%s'%(v,s) )

m.update()

# ======================== 3 Calculate expressions.
# ---expr_1.1---matrix__M_vs
M_vs = [ [0 for i in range(CNT_S)] for j in range(CNT_V) ]  # define the matrix__M_vs
for v in range(CNT_V):
for s in range(CNT_S):
M_vs[v][s] = dM_vs[v,s]
##print '\n----------M_vs:'
##print M_vs
##print '\n'

# --expr_1.2---dM_vk[v,k]
# --- Calaulate the M_vk rely on variable_matrix__M_vs.
M_vk = matrixMul(M_vs, M_sk)
##print '\n----------M_vk:'
##for element in M_vk:
##    print element
##print '\n'
dM_vk = {}
for v in range(CNT_V):
for k in range(CNT_K):
dM_vk[v,k] = M_vk[v][k]
##print '\n----------dM_vk:\n'
##for key, value in dM_vk.items():
##    print 'key:',key,' -- val:', value,'\n'
##print '---------------------------------\n'

# --expr_2--- placement relationship between VMs and Servers
f_ij_vs = {}
for i in range(CNT_V):
for j in range(CNT_V):
f_ij_vs[i,j]=quicksum(dM_vs[i,s]*dM_vs[j,s] for s in range(CNT_S))

# --expr_3--- placement relationship between VMs and Racks
f_ij_vk = {}
for i in range(CNT_V):
for j in range(CNT_V):
f_ij_vk[i,j]=quicksum(dM_vk[i,k]*dM_vk[j,k] for k in range(CNT_K))

# --expr_4--- communication cost C_ij between Vi and Vj
dC_ij = {}
for i,j in f_ij_vs:
dC_ij[i,j] = (1-f_ij_vk[i,j])*Cost_diffK + (f_ij_vk[i,j]-f_ij_vs[i,j])*Cost_diffS_sameK + f_ij_vs[i,j]*Cost_sameS

##print '\n----------dC_ij:\n'
##for key, value in dC_ij.items():
##    #print 'key:',key,' -- val:', value,'\n'
##    print 'key:',key,'\n'
##print '---------------------------------\n'

# ======================== 4 Add constriants:
# --- Constraint_1 ---- Completed Placement:
for i in range(CNT_V):
sumExpr_place_vs = quicksum( dM_vs[i,j] for j in range(CNT_S) )
m.addConstr( sumExpr_place_vs==1 )

# --- Constraint_2 ----  Resource Capacity:
# --- The Matrix_consum_rs: the consumption of Resource_s on Server_s
M_consum_rs = matrixMul( Q_rv, M_vs )
dM_consum_rs = {}
for r in range(CNT_R):
for s in range(CNT_S):
dM_consum_rs[r,s] = M_consum_rs[r][s]
##print '\n----------dM_consum_rs:\n'
##for key, value in dM_consum_rs.items():
##    print 'key:',key,' -- val:', value,'\n'
##print '---------------------------------\n'

for r in range(CNT_R):
for s in range(CNT_S):
m.addConstr( dM_consum_rs[r,s] <= Cap_rs[r][s] )

m.update()

# ======================== 5 set objective
sumCost = 0
for i in range(CNT_V):
for j in range(CNT_V):
sumCost = sumCost + T_ij[i][j]*dC_ij[i,j]
m.setObjective( sumCost )
m.update()

# ---- Compute optimal solution
m.optimize()

# ======================== 6 Print solution
Cost_OPT = 0
if m.status == GRB.status.OPTIMAL:
Cost_OPT = m.ObjVal
##        print '\n----------- Opt-sol-placement ------------:'
##        #print 'Opt.value =', m.ObjVal
##        for key,value in dM_vs.items():
##            print 'Pos:',key,' -- Yes/No:', value.getAttr("x")
##        print '---------------------------------:~\n'

else:
m.computeIIS()
m.write("model.ilp")

return Cost_OPT# Only return the OPT-value here, In fact, you can return all the solution via a dict.

2) C++ 嵌套调用 sample: ( 今晚时间有点紧，留到有空再贴这部分，2013-03-27)。（Ok，今天是三月最后一天，不太忙了终于，来把这个尾巴补上，2013-03-31。）

首先是用 Gurobi 解一个数学模型的通用逻辑的封装头文件：

/*   Deveploper: Davy_H
/*   2013-03-04
/*   @Dev: Embed Python_Gurobi into C++ project.
*/
#ifndef __OPT_MODULE_H__
#define __OPT_MODULE_H__

#include "Base_def.h"
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <set>
#include <fstream>
// -------------------------------- Embed Python_Gurobi into C++ project. _2013-03-04 ------------------------
#include <Python.h>
using namespace std;

class COPT_module
{
public:
COPT_module();
~COPT_module();

public:
// ----- 1. take over the argument.  // --- At last,  Prepare the inputing  str_Arg_for_model.
void Take_over_and_analyze_the_Inputing_data( _const_in_arg int nCNT_K, _const_in_arg int nCNT_S, _const_in_arg int nCNT_V, _const_in_arg int nCNT_R,
_const_in_arg Matrix & matT_ij, _const_in_arg Matrix & matCap_rs, _const_in_arg Matrix & matQ_rv, _const_in_arg Matrix & matM_sk );

// ----- 2. Begin to solve the model.
void Begin_to_solve_model( void );

// ----- 3. Get the optimal solution result.
int Get_optimal_value( void );

// ----- 4. Reset and Clear operations at the end of every solving.
void Reset_and_Clear_finally( void );

private:
int solve_model( void );
void set_opt_val( int nOpt_val );

private:
void reset_Opt_cost( void );
void Initialize_Python_gurobi_env( void );
void Finalize_Python_gurobi_env( void );

private:
int m_nOpt_cost;
string m_strInputingArg_into_model;
};

#endif

然后，只展示其中一个关键的方法，即 int solve_model(void); ，在这个方法里涉及到了关键的地方：在 c++ 里嵌套 python 脚本文件的方法，以及C++的数据结构
与 python数据结构 之间的转换。

int COPT_module::solve_model( void )
{
Initialize_Python_gurobi_env();      // Initialize the Python_Gurobi_environment firstly.

PyObject *pyMod = PyImport_ImportModule("vm_place"); 	// load the module, the "vm_place" is the name of the python_script file vm_place.py
if ( !pyMod )
{
cerr << " Err: pyMod == NULL " << endl;
return 0;
}
else  // if ok
{
//cout << " 1 ------- pyMod: " << pyMod << endl;
string strFunc_name = "SolveModel";	// The method name need to be specified.
PyObject *pyFunc = PyObject_GetAttrString(pyMod, strFunc_name.c_str());	    // load the function
if (pyFunc && PyCallable_Check(pyFunc))
{
PyObject *pyParams = PyTuple_New(1);    // create a python_tuple Parameter.
PyObject *pyInputArg = PyString_FromString(m_strInputingArg_into_model.c_str());  // Transform the C++ argument into Python_obj
PyTuple_SetItem(pyParams, 0, pyInputArg);

PyObject* pyRetObj = PyObject_CallObject(pyFunc, pyParams);	// ok, call the function in the model_Python_file
if ( ! pyRetObj  )
{
cerr << " Err: pyRetObj == NULL " << endl;
return 0;
}
else
{
//cout << " 2 ---------------- pyRetObj: " << pyRetObj << endl;
int nRet_cost_opt = 0;
nRet_cost_opt = PyInt_AsLong( pyRetObj );   // Translate the Python_integer to the C++ long int.

if ( nRet_cost_opt > 0 )    // Return the optimal value.
return nRet_cost_opt;
}
}
else
{
cerr << " Err: pyParams == NULL " << endl;
return 0;
}
}

return 0;
}

简单解释：这里我只是返回了一个整型数，所以无需对 PyObject_CallObject 的返回值进行解析，只需简单地用 PyInt_AsLong()
函数把Python int 转换为C++ 的 long int 即可。否则，如果从
Python写的数学模型程序 中返回的是一个 python dict 的数据结构，那么就需要用其他的相关 api 来解析了，对更复杂结构的转化的例子有空再贴。

-- Davy Hwang

2013-03-31

// ============================ //

等了一年多，现在是2015年01月23号，今天又用到这一处的数据转换，我索性抽出一点时间，把这最后一个尾巴补上吧。这里我举几个例子，来演示在 C++ 程序中对 从 Python写的数学模型程序 中返回数据结构的解析方法与过程。

例子1：如果从用Python写的数学模型程序中返回的只是一个整型数或者浮点数，那么在c++这头是如何接收并转换为double类型的。

1.1）python 程序的返回值部分代码：

def Func_solve( Input_Args ):
# ... do something
return nCost

1.2）C++这边接收并转换数据类型：

double COPT_module::solve_math_model( void )
{
Initialize_Python_gurobi_env();
PyObject *pyMod = PyImport_ImportModule("Name_of_python_model");	// load the name of the math-model file.

if ( pyMod == NULL )
{
cout << " 1 -------err--------- pyMod == NULL " << endl;
return 0;
}

else  // if ok
{
string strFunc_name = "Func_solve";	// This is the Function_name_in_the_math_model
PyObject *pyFunc = PyObject_GetAttrString(pyMod, strFunc_name.c_str());	// load the function
if ( ! pyFunc )
{
cout << " 2 -------err--------- pyFunc == NULL " << endl;
return 0;
}

int nCallable_check = PyCallable_Check( pyFunc );
if (  nCallable_check <= 0 )
{
cout << " 3 -------err--------- pyFunc is not Callable. " << endl;
return 0;
}
else
{
PyObject *pyParams = PyTuple_New(1);
if ( !pyParams )
{
cout << " 4 -------err--------- pyParams == NULL " << endl;
return 0;
}

// m_strInputingArg_into_trans_model is the Inputing "Input_Args" to the math-model.
PyObject *pyInputArg = PyString_FromString( m_strInputingArg_into_trans_model.c_str() );
PyTuple_SetItem(pyParams, 0, pyInputArg);

PyObject* pyRetObj = PyObject_CallObject(pyFunc, pyParams);	// ok, call the function
if ( !pyRetObj )
{
cout << " 5 -------err--------- pyRetObj == NULL " << endl;
return 0;
}
else
{
// ===================== Get the returned_value
double dRet_cost_returned = PyFloat_AsDouble( pyRetObj );
return dRet_cost_returned;
}
}
}

return 0;
}

例子2：如果从用Python写的数学模型程序中返回的是一个字典结构，那么在c++这头是如何接收并解析字典中的返回值的。

2.1）python 程序的返回 dict 的部分代码：

def Func_solve( Input_Args ):
# ... do something
dict_Cost_returned = {}	### Suppose that there will be two elements in this returned dict.
dict_Cost_returned[0] = nCost_1st
dict_Cost_returned[1] = nCost_2nd
return dict_Cost_returned

2.2）C++这边接收并转换 dict 数据类型（还是在上面那段代码的基础上，只需要修改最关键的else{}那一段内容）：
//...

PyObject* pyRetObj = PyObject_CallObject(pyFunc, pyParams);	// ok, call the function
if ( !pyRetObj )
{
cout << " 5 -------err--------- pyRetObj == NULL " << endl;
return;
}
else
{
double dRet_cost_1 = 0;
double dRet_cost_2 = 0;

// ==================== 0. Get the keys of the returned dict.
PyObject* pKeys = PyDict_Keys( pyRetObj );

// -------------------- 1. get and Save the 1st returned value.
PyObject* pCol_0 = PyList_GetItem( pKeys, 0 );
PyObject* pValues_cost_0 = PyDict_GetItem( pyRetObj, pCol_0 );
dRet_cost_1 = PyFloat_AsDouble( pValues_cost_0 );
// save dRet_cost_1

// -------------------- 2  get and  Save the 2nd returned value.
PyObject* pCol_1 = PyList_GetItem( pKeys, 1 );
PyObject* pValues_cost_1= PyDict_GetItem( pyRetObj, pCol_1 );
dRet_cost_2 = PyFloat_AsDouble( pValues_cost_1 );
// save dRet_cost_2
}

//...

例子3：如果从用Python写的数学模型程序中返回的是一个字典结构，而且字典里有两个元素：第一个是个整数，第二个是含有多个整数的list.

3.1）python返回程序省略（请各种发挥想象）。
3.2）C++这边接收并转换解析 dict, list 数据类型（还是在第一段代码的基础上，只需要修改最关键的else{}那一段内容）：

// ...
PyObject* pyRetObj = PyObject_CallObject(pyFunc, pyParams);	// ok, call the function
if ( !pyRetObj )
{
cout << " 5 -------err--------- pyRetObj == NULL " << endl;
return 0;
}
else
{
// ===================== Get the keys of the returned_dict.
PyObject* pKeys = PyDict_Keys( pyRetObj );

// -------------------- 1. get and Save the 1st returned double-element.
PyObject* pCol_0 = PyList_GetItem( pKeys, 0 );
double dRet_cost_1st_element = 0;
PyObject* pValues_cost_opt = PyDict_GetItem( pyRetObj, pCol_0 );
dRet_cost_1st_element = PyFloat_AsDouble( pValues_cost_opt );

// -------------------- 2  get and  Save the second returned element: a list.
PyObject* pCol_1 = PyList_GetItem( pKeys, 1 );
PyObject* pList_returned_2nd_element = PyDict_GetItem( pyRetObj, pCol_1 );
Py_ssize_t nlen_List_returned = PyObject_Length( pList_returned_2nd_element );

for ( int i = 0; i < nlen_List_returned; ++i )
{
PyObject* pVal_in_list = PyList_GetItem( pList_returned_2nd_element, i );
int nVal_i = PyInt_AsLong( pVal_in_list );
m_vector_values_in_the_returned_list.push_back( nVal_i );
}
}
// ...

经过以上三个例子，相信所有的情况都可以套用它们来解决。至此，这篇blog就完整了，希望某一天对看到它的某位读者有参考意义。

Davy_H
2015-01-23