
import numpy as np
import gurobipy as gp

c = np.array([4,5])
h = np.array([2,-7,5])
F = np.concatenate((np.array([[3,4]]),np.zeros((3,2))),axis=0)
G = np.concatenate((np.array([[2,-2,3]]),np.eye(3)),axis=0)
d = np.array([10,2,2,2])

x_star=np.array([3,3])


def solve_DSP(d,h,F,G,xstar,silent=False):
    # Model
    m = gp.Model("DSP")
    m.setParam(gp.GRB.Param.InfUnbdInfo, 1)
    if silent:
        m.setParam(gp.GRB.Param.OutputFlag, 0)
    
    # Create variables
    u = m.addMVar(shape=G.shape[0], vtype=gp.GRB.CONTINUOUS, name="u")

    # Set objective
    m.setObjective(u @ (d-F @ xstar), gp.GRB.MINIMIZE)

    # Add constraints
    m.addConstr(G.T @ u>= h, name="c")
    if not silent: 
        m.update()
        m.display()
    # Optimize model
    m.optimize()
    return m, u


dsp_model, dsp_sol = solve_DSP(d,h,F,G,x_star)

if dsp_model.status == gp.GRB.status.INFEASIBLE:
    print("Instance infeaasible")
elif dsp_model.status == gp.GRB.status.OPTIMAL:
    print("extreme point:", dsp_sol.X)
elif dsp_model.status == gp.GRB.status.UNBOUNDED:
    print("extreme ray:",dsp_sol.UnbdRay)

Set parameter Username
Academic license - for non-commercial use only - expires 2024-02-26
Set parameter InfUnbdInfo to value 1
Minimize
  -11.0 u[0] + 2.0 u[1] + 2.0 u[2] + 2.0 u[3]
Subject To
  c[0]: 2.0 u[0] + u[1] >= 2
  c[1]: -2.0 u[0] + u[2] >= -7
  c[2]: 3.0 u[0] + u[3] >= 5
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (linux64)

CPU model: Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz, instruction set [SSE2|AVX]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 3 rows, 4 columns and 6 nonzeros
Model fingerprint: 0x91a40911
Coefficient statistics:
  Matrix range     [1e+00, 3e+00]
  Objective range  [2e+00, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [2e+00, 7e+00]
Presolve time: 0.01s
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -1.1000000e+31   1.000000e+30   1.100000e+01      0s

Solved in 2 iterations and 0.01 seconds (0.00 work units)
Unbounded model
extreme ray: [0.5 0.  1.  0. ]


dsp_sol.UnbdRay @ (d-F @ x_star)

-3.5


def solve_REF(c,d,F,extreme_rays=[], extreme_points=[], silent=False):
     # Model
    m = gp.Model("DSP")
    m.setParam(gp.GRB.Param.DualReductions, 0)
    if silent:
        m.setParam(gp.GRB.Param.OutputFlag, 0)

    # Create variables
    x = m.addMVar(shape=F.shape[1], ub=3, vtype=gp.GRB.CONTINUOUS, name="x")
    eta = m.addMVar(shape=1, ub=100, vtype=gp.GRB.CONTINUOUS, name="eta")

    # Set objective
    m.setObjective(c @ x + eta, gp.GRB.MAXIMIZE)

    # Add constraints
    for ray in extreme_rays:
        m.addConstr(ray @ (d - F @ x)>=0, name="feasibility cut")
    for point in extreme_points:
        m.addConstr(point @ (d - F @ x)>=eta, name="optimality cut")
    if not silent:
        m.update()
        m.display()
    # Optimize model
    m.optimize()
    return m,x,eta


ref_model, ref_sol_x, ref_sol_eta = solve_REF(c,d,F, [dsp_sol.UnbdRay])

Set parameter DualReductions to value 0
Maximize
  4.0 x[0] + 5.0 x[1] + eta[0]
Subject To
  feasibility cut: -1.5 x[0] + -2.0 x[1] >= -7
Bounds
  0 <= x[0] <= 3
  0 <= x[1] <= 3
  0 <= eta[0] <= 100
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (linux64)

CPU model: Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz, instruction set [SSE2|AVX]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 1 rows, 3 columns and 2 nonzeros
Model fingerprint: 0x8bff916f
Coefficient statistics:
  Matrix range     [2e+00, 2e+00]
  Objective range  [1e+00, 5e+00]
  Bounds range     [3e+00, 1e+02]
  RHS range        [7e+00, 7e+00]
Presolve time: 0.01s
Presolved: 1 rows, 3 columns, 2 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.1866667e+02   8.333333e-01   0.000000e+00      0s
       1    1.1825000e+02   0.000000e+00   0.000000e+00      0s

Solved in 1 iterations and 0.01 seconds (0.00 work units)
Optimal objective  1.182500000e+02


extreme_rays=[]
extreme_points=[]
while True:
    # We solve the Restricted Benders reformulation (RBR) or restricted extended formulation:
    ref_model, ref_sol_x, ref_sol_eta = solve_REF(c,d,F,extreme_rays, extreme_points,silent=True)
    if ref_model.status == gp.GRB.status.INFEASIBLE:
        print("Instance infeaasible")
        break
    elif ref_model.status == gp.GRB.status.UNBOUNDED:
        print("Ref unbounded:",ref_sol_x.X,ref_sol_eta.X)        
        break
    elif ref_model.status == gp.GRB.status.OPTIMAL:
        print("REF solution",ref_sol_x.X,ref_sol_eta.X)
    else:
        print('Optimization was stopped with status %d' % m.status)
        break
    # We solve the Dual Subproblem
    dsp_model, dsp_sol = solve_DSP(d,h,F,G,ref_sol_x.X, silent=True)
    if dsp_model.status == gp.GRB.status.UNBOUNDED:
        print("Found extreme ray:",dsp_sol.UnbdRay)
        extreme_rays=extreme_rays+[dsp_sol.UnbdRay]
    elif dsp_model.status == gp.GRB.status.OPTIMAL:
        print("Found extreme point:", dsp_sol.X)
        if dsp_model.objVal < ref_sol_eta.X:
            extreme_points=extreme_points+[dsp_sol.X]
        elif dsp_model.objVal == ref_sol_eta.X:
            print("Problem Solved: ",
                  f"z={ref_model.objVal}, eta={ref_sol_eta.X},"
                  f"x={ref_sol_x.X}, y={[c.Pi for c in dsp_model.getConstrs()]}")
            break
    elif dsp_model.status == gp.GRB.status.INFEASIBLE:
        print("DSP Instance infeaasible")
        break

Set parameter DualReductions to value 0
REF solution [3. 3.] [100.]
Set parameter InfUnbdInfo to value 1
Found extreme ray: [0.5 0.  1.  0. ]
Set parameter DualReductions to value 0
REF solution [3.   1.25] [100.]
Set parameter InfUnbdInfo to value 1
Found extreme point: [3.5 0.  0.  0. ]
Set parameter DualReductions to value 0
REF solution [0. 0.] [35.]
Set parameter InfUnbdInfo to value 1
Found extreme point: [0. 2. 0. 5.]
Set parameter DualReductions to value 0
REF solution [2. 0.] [14.]
Set parameter InfUnbdInfo to value 1
Found extreme point: [1.66666667 0.         0.         0.        ]
Set parameter DualReductions to value 0
REF solution [0.53333333 0.        ] [14.]
Set parameter InfUnbdInfo to value 1
Found extreme point: [1. 0. 0. 2.]
Set parameter DualReductions to value 0
REF solution [1.33333333 0.        ] [10.]
Set parameter InfUnbdInfo to value 1
Found extreme point: [1. 0. 0. 2.]
Problem Solved:  z=15.333333333333334, eta=[10.],x=[1.33333333 0.        ], y=[-8.881784197001252e-16, 0.0, 2.0]

Exercises on Benders Algorithm¶

Task 1¶

Task 2¶

Solution¶