# 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