# Experiments for the small inventory example
# Frans de ruiter, 2017 creative commons.
using JuMP, Mosek
# Input parameters
const T = 2 # Number of periods
const Γ = 5 # Size of uncertainty set
const h = 1 # Holding costs
const b = 2 # Backlog costs
const I_0 = 5 # Initial inventory costs
const dnom = 5*ones(T) # nominal demand
const qmax = 3 # capacity in second period
helpmatr = [i<=j ? 1 : 0 for i=1:T, j=1:T] # upper triangular with ones for nonanticipativity rule
# Declare model
m = Model(solver=MosekSolver())
# Variables
@variable(m, obj)
@variable(m, c_const[1:T])
@variable(m, c_coeff[1:T,1:3,1:T])
@variable(m, q_const[1:T])
@variable(m, q_coeff[1:T,1:3,1:T])
# Nonanticipativity: require entries of q_coeff[t,:,:] to be zero if they correspond to demand from current or future periods
@constraint(m, [i=1:T,j=1:3, s=i:T], q_coeff[i,j,s] == 0)
# !!!! OUTCOMMENT constraints below for pure linear decision rules
# @constraint(m, q_coeff[:,2:3,:] .== 0)
# @constraint(m, c_coeff[:,2:3,:] .== 0)
# !!!!
# !!!! OUTCOMMENT constraints below for static robust solution for c_1
# @constraint(m, c_coeff[:,1:3,:] .== 0)
# !!!!
# @constraint(m, q_const[1]==0)
# Function to add the support function used to reformulate the robust counterpart using the lifted set V
function support(w,wabs,wsq)
ϕ = @variable(m, lowerbound = 0)
η = @variable(m, [1:T])
t = @variable(m, [1:T], lowerbound = 0)
@constraint(m, [i=1:T], norm( [wabs[i]+η[i], ϕ - wsq[i]-t[i]] ) <= ϕ - wsq[i] + t[i])
@constraint(m, [i=1:T], w[i] <= η[i])
@constraint(m, [i=1:T], -w[i] <= η[i])
@constraint(m, [i=1:T], wsq[i] <= ϕ)
return (Γ^2)*ϕ + sum(t)
# return Γ*norm(w)
end
# Objective
@objective(m, Min, obj)
# Constraint for objective value
@constraint(m, sum(c_const) + support( sum(c_coeff[s,1,:] for s=1:T), sum(c_coeff[s,2,:] for s=1:T), sum(c_coeff[s,3,:] for s=1:T) ) <= obj)
# Constraints for holding costs
@constraint(m, [s=1:T], -(1/h)*c_const[s] + I_0 + sum(q_const[1:s]) - sum(dnom[1:s]) + support( -(1/h)*c_coeff[s,1,:] + sum(q_coeff[τ,1,:] for τ=1:s)-helpmatr[:,s], -(1/h)*c_coeff[s,2,:] + sum(q_coeff[τ,2,:] for τ=1:s), -(1/h)*c_coeff[s,3,:] + sum(q_coeff[τ,3,:] for τ=1:s)) <= 0)
# Constraints for backlog costs
@constraint(m, [s=1:T], -(1/b)*c_const[s] - I_0 - sum(q_const[1:s]) + sum(dnom[1:s]) + support( -(1/b)*c_coeff[s,1,:] - sum(q_coeff[τ,1,:] for τ=1:s) + helpmatr[1:T,s], -(1/b)*c_coeff[s,2,:] - sum(q_coeff[τ,2,:] for τ=1:s), -(1/b)*c_coeff[s,3,:] - sum(q_coeff[τ,3,:] for τ=1:s)) <= 0)
# Capacity constraint
@constraint(m, q_const[2] + support(q_coeff[2,1,:],q_coeff[2,2,:],q_coeff[2,3,:]) <= qmax)
# Nonnegativity constraints
@constraint(m, [s=1:T], -q_const[s] + support(-q_coeff[s,1,:],-q_coeff[s,2,:],-q_coeff[s,3,:]) <= 0)
#
solve(m)
# Solution:
q0 = getvalue(q_const)
q1 = getvalue(q_coeff)
c0 = getvalue(c_const)
c1 = getvalue(c_coeff)
WCobj = getvalue(obj)
# ! Note that the linear decision rule is in terms of ζ. Hence, the linear decision rule
#
# q(ζ) = a + bζ + c|ζ|,
#
# is equivalent to (using the transformation d = ζ + 5)
# q(d) = (a-5b) + bd + c|d-5|
#
#
# End of File