用python做整数规划

import pyomo.environ as pyo

def build_model(v, a=0):
    model = pyo.ConcreteModel()

    # creat indices and sets
    model.I = pyo.Set(initialize=range(1,51))
    model.J = pyo.Set(initialize=list(range(1,201)))

    # #input data
    # mu = 5
    # sigma = 10
    # def c_init(model, i, j):
    #     c_data = np.random.randn(1) * sigma +mu
    # model.f = pyo.Param(model.I, initialize = 50)
    # model.c = pyo.Param(model.I, model.J, initialize = c_init)

    # input data
    model.f = pyo.Param(model.I, initialize = 50)
    model.c = pyo.Param(model.I, model.J, initialize=v)

    # decision varibles
    model.x = pyo.Var(model.I, model.J, within=pyo.Binary)  # give a initial value is crucial
    model.y = pyo.Var(model.I, within=pyo.Binary)

    # objective function
    def Obj_rule(m):
        return pyo.summation(m.f, m.y) + sum(sum(m.c[i,j] * m.x[i,j] for j in m.J) for i in m.I)

    model.Obj = pyo.Objective(rule=Obj_rule)

    # Constraints
    def must_be_served(m, j):
        return sum(m.x[i,j] for i in m.I) == 1

    model.must_be_served_constraints = pyo.Constraint(model.J, rule=must_be_served)

    if a == 0:
        def facility_open(m, i):
            return  sum(m.x[i,j] for j in m.J) <= 100*m.y[i]
        model.facility_open_constraint = pyo.Constraint(model.I,rule=facility_open)
    else:
        def facility_open(m, i,j):
            return m.x[i,j]<= m.y[i]
        model.facility_open_constraint = pyo.Constraint(model.I, model.J, rule=facility_open)

import random
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

import time 

solver = SolverFactory('coptampl')

import Lab5_homework_model1_Q4 as Lab5_homework_model1

def v_init():
    v = dict()
    for i in list(range(1,51)):
        for j in list(range(1,201)):
            #v[i,j] = np.random.randn(1)* 5 + 10
            v[i,j] = random.random()*5 + 10
    return v

for itr in range(5):
    v = v_init()
    
    model1 = Lab5_homework_model1.build_model(v,a=0)
    
    model2 = Lab5_homework_model1.build_model(v,a=1)
    
    time_st1 = time.time()
    results1 = solver.solve(model1, tee = True)
    print(results1)
    time_end1 = time.time()
    print("The solver returned a status of:" + str(results1.solver.status))

    results2 = solver.solve(model2, tee = True)
    print(results2)
    time_end2 = time.time()
    print("The solver returned a status of:" + str(results2.solver.status))

    with open('Q4.txt','a') as f:
        f.write('ObjectVal of model1 for ' + str(itr+1) + ' times result:' +str(pyo.value(model1.Obj)) + '  used time(sec): ' + str(time_end1-time_st1))
        f.write('\r')
        f.write('ObjectVal of model2 for ' + str(itr+1) + ' times result:' +str(pyo.value(model2.Obj)) + '  used time(sec): ' + str(time_end2-time_end1))
        f.write('\r')