import numpy as np
import gurobipy as gp

# define Gurobi model
m = gp.Model("Unit Commitment")

n = 3
ubs = [100, 90, 80]
demands = [185, 150, 150]
lb = 90

# define variables
x = m.addMVar(n, ub=ubs, name="x")

# add constraints and names
m.addConstr(sum(x[0], x[1]) >= demands[0], name="demand_0")
m.addConstr(sum(x[0], x[2]) >= demands[1], name="demand_1")
m.addConstr(sum(x[1], x[2]) >= demands[2], name="demand_2")
m.addConstr(x[0] >= lb, name="lb")

# define objective function
c = np.array([10, 8, 9])
m.setObjective(c@x, gp.GRB.MINIMIZE)

# solve
m.optimize()

# check status
if m.status == gp.GRB.INFEASIBLE:
    print("infeasible")
elif m.status == gp.GRB.OPTIMAL:
    print("feasible")

    x_sol = x.X
    y_sol = np.array([
        m.getAttr("Pi", [m.getConstrByName("demand_{}".format(i))])[0]
        for i in range(n)
    ])
    val = m.objVal
print("x={} y={} obj={}\n".format(x_sol, y_sol, val))

# gather existing constraints
constrs_list = []
for cs in m.getConstrs():
    sense = cs.sense
    name = cs.ConstrName
    if "demand_" in name:
        constrs_list.append(cs)
print(constrs_list)
col = gp.Column(np.ones(len(constrs_list)), constrs_list)

# add new variable
x3 = m.addVar(obj=8, name='x[3]', column=col)
m.update()
x = m.getVars()

# add new upper bound constraints
m.addConstr(x[3] <= 20)

# modify RHS
m.setAttr("RHS", m.getConstrByName("lb"), 50)

# re-solve and check
m.optimize()
if m.status == gp.GRB.INFEASIBLE:
    print("infeasible")
elif m.status == gp.GRB.OPTIMAL:
    print("feasible")

    x_sol = [x[i].X for i in range(len(x))]
    y_sol = np.array([
        m.getAttr("Pi", [m.getConstrByName("demand_{}".format(i))])[0]
        for i in range(n)
    ])
    val = m.objVal

print("x={} y={} obj={}\n".format(x_sol, y_sol, val))

m.setParam('OutputFlag', 0)
m.setParam(gp.GRB.Param.Threads, 4)
m.setParam(gp.GRB.Param.TimeLimit, 60*60) # 1h time limit
m.setParam(gp.GRB.Param.Crossover, 0)
m.update()

m.optimize()
if m.status == gp.GRB.INFEASIBLE:
    print("infeasible")
elif m.status == gp.GRB.OPTIMAL:
    print("feasible")

    x_sol = [x[i].X for i in range(len(x))]
    y_sol = np.array([
        m.getAttr("Pi", [m.getConstrByName("demand_{}".format(i))])[0]
        for i in range(n)
    ])
    val = m.objVal

print("x={} y={} obj={}\n".format(x_sol, y_sol, val))