docplex pyscipopt 精确求解例子

尝试过了群智能算法和邻域搜索，精确求解得来个栗子才平衡，是不?

假设鸡肉，白菜两种食物的价格 和 某些维生素含量， 以及最小摄入量的关系如下表：
请问两种食物 各食用多少，既满足需要，费用又最省？

	鸡肉	白菜	最低摄入量
维生素I	0.1	0.15	1
维生素J	1.7	0.75	7.5
维生素K	1.1	1.3	10
单价	2	1.5

最优化问题， 一般有三要素：
1 决策变量  通常就是要求解的未知量x
2 目标函数  通常是 要优化(最大或最小)的那个目标的数学表达式  是决策变量的函数 f(x)
3 约束条件  对决策变量的限制条件 即 x的允许取值范围  称为可行域
 

所以看得到代码中，求解的步骤基本也是：
增加决策变量
增加约束
设置目标函数
求解

线性规划的标准形(矩阵形式)长这样：

约束条件为等式 且b>=0
决策变量x非负

如果约束条件不为等式怎么办？ 用松驰变量 或 剩余变量 转化为等式。

from collections import namedtuple
import operator
from pyscipopt import quicksum
from pyscipopt import Model as ScipModel
from docplex.mp.model import Model as CplexModel


# min 2 * x1 + 1.5 * x2
# subject to
#     0.1 * x1 + 0.15 * x2 >= 1.0
#     1.7 * x1 + 0.75 * x2 >= 7.5
#     1.1 * x1 + 1.30 * x2 >=10.0
#     x1 >= 0
#     x2 >= 0



QuestionData = namedtuple("QuestionData", ("unit_price", "low_b", "var_types",
                                           "const_A", "operator_symbol", "const_b"))


class ExactSolution:
    def __init__(self):
        pass

    def __call__(self):
        self.add_decision_var()
        self.add_constraint()
        self.set_obj()
        # start solve
        self.optimize()
        return None

    def add_decision_var(self):
        raise NotImplementedError

    def add_constraint(self):
        raise NotImplementedError

    def set_obj(self):
        raise NotImplementedError

    def optimize(self):
        raise NotImplementedError

    def get_result(self):
        raise NotImplementedError


class ScipSolution(ExactSolution):
    def __init__(self, model_name, question_data):
        self.model = ScipModel(model_name)
        self.question_data = question_data
        super().__init__()

    def add_decision_var(self):
        self.x_var = {}
        for ix, (var_lb, var_type) in enumerate(zip(self.question_data.low_b, self.question_data.var_types)):
            self.x_var[ix] = self.model.addVar(vtype=var_type, lb=var_lb)

    def add_constraint(self):
        for const_a_group, op, b in zip(
                self.question_data.const_A, self.question_data.operator_symbol, self.question_data.const_b):
            self.model.addCons(op(quicksum(const_a * self.x_var[ix]
                               for const_a, ix in zip(const_a_group, list(range(len(self.x_var))))), b))

    def set_obj(self):
        self.model.setObjective(quicksum(u_price * self.x_var[ix]
                                         for u_price, ix in zip(self.question_data.unit_price, range(len(self.x_var)))),
                                sense="minimize")

    def optimize(self):
        self.model.optimize()

    def get_result(self):
        obj = self.model.getObjVal()
        x_value = [self.model.getVal(self.x_var[ix]) for ix in range(len(self.x_var))]
        print("\nScipopt: the result cost is %s and the quantity of vector x is %s " % (obj, x_value))
        return obj, x_value


class CplexSolution(ExactSolution):
    def __init__(self, model_name, question_data):
        self.model = CplexModel(model_name)
        self.question_data = question_data
        super().__init__()

    def add_decision_var(self):
        self.x_var = {}
        for ix, (var_lb, var_type) in enumerate(zip(self.question_data.low_b, self.question_data.var_types)):
            if var_type == "C":
                self.x_var[ix] = self.model.continuous_var(lb=var_lb)
            elif var_type == "I":
                self.x_var[ix] = self.model.integer_var(lb=var_lb)
            elif var_type == "B":
                self.x_var[ix] = self.model.binary_var()
            else:
                raise ValueError("the var type %s is wrong!" % var_type)

    def add_constraint(self):
        self.model.add_constraints(op(sum(const_a * self.x_var[ix]
                                      for const_a, ix in zip(const_a_group, list(range(len(self.x_var))))), b)
                                   for const_a_group, op, b in zip(
            self.question_data.const_A, self.question_data.operator_symbol, self.question_data.const_b))


    def set_obj(self):
        # self.model.minimize(self.model.sum(self.e1_var[s_id] + self.e2_var[s_id] for s_id in self.sample_no))
        self.model.minimize(sum(u_price * self.x_var[ix]
                            for u_price, ix in zip(self.question_data.unit_price, range(len(self.x_var)))))


    def optimize(self):
        self.solution = self.model.solve(log_output=False)
        if not self.solution:
            raise RuntimeError("can't get solution!")

    def get_result(self):
        obj = self.solution.objective_value
        x_value = [self.x_var[ix].solution_value for ix in range(len(self.x_var))]
        print("\nCplex: the result cost is %s and the quantity of vector x is %s " % (obj, x_value))
        return obj, x_value


def optimiz(question_data, solver_type='cplex', model_name='food_sample'):
    if solver_type == 'cplex':
        sol = CplexSolution(model_name, question_data)
        sol()
    elif solver_type == 'scip':
        sol = ScipSolution(model_name, question_data)
        sol()
    else:
        raise ValueError("the solver_type %s is not allowed!" % solver_type)
    return sol


def get_question_data():
    # x1, x2
    unit_price = [2.0, 1.5]
    low_b = [0.0, 0.0]
    var_types = 'CC'
    const_A = [[0.1, 0.15], [1.7, 0.75], [1.1, 1.3]]
    operator_symbol = [operator.ge, operator.ge, operator.ge]
    const_b = [1.0, 7.5, 10.0]
    return QuestionData(unit_price, low_b, var_types, const_A, operator_symbol, const_b)


def run(solver_type='cplex'):
    question_data = get_question_data()
    sol = optimiz(question_data, solver_type=solver_type)
    obj, x_value = sol.get_result()


if __name__ == '__main__':
    run(solver_type='cplex')     # "scip" or "cplex"

运行结果：

Cplex: the result cost is 12.725631768953068 and the quantity of vector x is [1.6245487364620947, 6.317689530685919]

Scipopt: the result cost is 12.725631768953068 and the quantity of vector x is [1.6245487364620947, 6.317689530685919]