5 Master Problem

1, Functions in FuncMas.jl

function setModelMas(numQ, vecP, numSub, vecSenseP, dualPen, gurobi_env, whiSolver)
    vecDualGuessPi = 100 * ones(Float64, numQ)
    vecDualGuessKappa = 0 * ones(Float64, numSub)
    ## X1: initial matrix of extreme points
    if whiSolver == 1
        modMas =  Model(solver = GurobiSolver(OutputFlag = 0, gurobi_env))
    elseif whiSolver == 2
        modMas =  Model(solver = CplexSolver(CPX_PARAM_SCRIND=0))
    else
        modMas =  Model(solver = GLPKSolverLP())
    end
    K = 1
    # In this case we do not use a starting set of extreme points.
    # we just use a dummy starting column
    @variable(modMas, vecLambda[1:K] >= 0 )
    @variable(modMas, 0 <= vecMuMinus[1:numQ] <= dualPen)
    @variable(modMas, - dualPen <= vecMuPlus[1:numQ] <= 0)
    @variable(modMas, 0 <= vecMuMinusConv[1:numSub] <= dualPen)
    @variable(modMas, - dualPen <= vecMuPlusConv[1:numSub] <= 0)
    # Remember to consider if we need to maximize or minimize
    @objective(
        modMas, Max, sum(- 1000 * vecLambda[j] for j = 1:K) +
        sum(vecDualGuessPi[i] * (vecMuMinus[i] + vecMuPlus[i]) for i = 1:numQ) +
        sum(vecDualGuessKappa[k] * (vecMuMinusConv[k] + vecMuPlusConv[k]) for k=1:numSub)
        )
    # @constraint(modMas, vecConsRef[i = 1: numQ], sum(vecP[i] * vecLambda[j] for j = 1:K) == vecP[i])
    JuMPConstraintRef = JuMP.ConstraintRef{Model, JuMP.GenericRangeConstraint{JuMP.GenericAffExpr{Float64,Variable}}}
    vecConsRef = Vector{JuMPConstraintRef}(undef, numQ)
    for i = 1:numQ
        if vecSenseP[i] == leq
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] <= vecP[i])
        elseif vecSenseP[i] == geq
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] >= vecP[i])
        else
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] == vecP[i])
        end
    end
    @constraint(modMas, vecConsConvex[k = 1: numSub], sum(vecLambda[j] for j = 1: K) +
        vecMuMinusConv[k] + vecMuPlusConv[k] == 1)
    return (modMas, vecConsRef, vecConsConvex, vecLambda, vecMuMinus, vecMuPlus, vecMuMinusConv, vecMuPlusConv)
end

function solveMas(modMas, vecConsRef, vecConsConvex)
    status = solve(modMas)
    if status != :Optimal
        throw("Error: Non-optimal modMas-problem status")
    end
    vecPi = getdual(vecConsRef)
    vecKappa = getdual(vecConsConvex)
    obj_master = getobjectivevalue(modMas)
    ## ensure that vecPi and vecKappa are row vectors
    vecPi = reshape(vecPi, 1, length(vecPi))
    vecKappa = reshape(vecKappa, 1, length(vecKappa))
    return (vecPi, vecKappa, obj_master)
end