Linear programming with Jump

Specific Domains Optimization (Mathematical)
1

I am trying to write a linear programming problem problem in JuMP but I am getting error messages. I have given my code below, and I believe the issue is in definition of the constraint. Please can someone help to fix the code?

Apologies the format makes it difficult to read the code, but I cannot see an option to send the code as an attachment.

using JuMP
using Cbc

Preparing an optimization model

myModel = Model(with_optimizer(Cbc.Optimizer))

prices = [99.74 91.22 98.71 103.75 97.15]
cashFlows = [4 5 2.5 5 4 4 5 2.5 5 4 4 5 2.5 5 4 4 5 2.5 5 4 4 5 102.5 5 4 4 5 0 105 104 4 105 0 0 0 104 0 0 0 0]
Liab_CFs = [5 7 7 6 8 7 20 0]‘*1000
nt=size(cashFlows,1)
nb=size(cashFlows,2)
Rates = [0.01 0.015 0.017 0.019 0.02 0.025 0.027 0.029]’

Disc= [0.99009901 0.970661749 0.950686094 0.927477246 0.90573081 0.862296866 0.829863942 0.795567442]

#Number of bonds available
nBonds = [10 100 20 30 5]’

Declaring variables

@variable(myModel, 0<= x <=nBonds)

Setting the objective

The objective is to determine the cheapest set of bonds to cover liability cash flows subject to the constraint below.

@objective(myModel,Min,prices*x)

Add constraint

The constraint is determined by cash flows shortfall calculated as net cash flows (bonds cash flows - liability cash flows)

compounded at forward rate (B), and discounted to time zero. The maximum shortfall should be less than or equal to 0.05*53.844.

The max is linearized by setting each element in the matrix to be less than or equal to 0.05*53.844

B=[1.01 1.020024752 1.02101183 1.02502363 1.024009823 1.050370059 1.039082218 1.043109481]’
A=-(cashFlowsx-Liab_CFs)
d=cumprod(B[1:end-1])
d=[1 d[:]]';
M=tril(d./transpose(d))
@constraint(myModel, constraint1, ((MA).Disc)<=0.0553.844)

M*A is doing the following operation

M*A=[A[1];A[1]*B[1]+A[2];(A[1]*B[1]+A[2])*B[2]+A[3];(A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4];…

((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5]

(((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6];…

((((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6])*B[6]+A[7];…

(((((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6])*B[6]+A[7])*B[7]+A[8]]

Printing the prepared optimization model

print(myModel)

Solving the optimization problem

optimize!(myModel)

Printing the optimal solutions obtained

println(“Optimal Solutions:”)
println("x = ", JuMP.value(x))

2

Please see Please read: make it easier to help you. It describes how you format your code.

The likely issue is that you need:

@variable(myModel, 0<= x[b = 1:length(nBonds)] <= nBonds[b])

Here is the relevant documentation: Variables · JuMP

1 Like
3

Thank you for your help. I have made the change you suggested, but I am still getting error messages (see below). I have been through the documentation but I have not been able to fix this :

ERROR MESSAGE
The objective function GenericAffExpr{Float64,VariableRef}[99.74 x[1] + 91.22 x[2] + 98.71 x[3] + 103.75 x[4] + 97.15 x[5]] is not supported by JuMP.

Stacktrace:
[1] error(::String) at .\error.jl:33
[2] set_objective(::Model, ::MathOptInterface.OptimizationSense, ::Array{GenericAffExpr{Float64,VariableRef},1}) at C:\Users\a.julia\packages\JuMP\MsUSY\src\objective.jl:101
[3] top-level scope at C:\Users\a.julia\packages\JuMP\MsUSY\src\macros.jl:979
[4] top-level scope at In[6]:21

CODE

using JuMP
using Cbc

# Preparing an optimization model
myModel = Model(with_optimizer(Cbc.Optimizer))

prices = [99.74 91.22 98.71 103.75 97.15]
cashFlows = [4 5 2.5 5 4 4 5 2.5 5 4 4 5 2.5 5 4 4 5 2.5 5 4 4 5 102.5 5 4 4 5 0 105 104 4 105 0 0 0 104 0 0 0 0]
Liab_CFs = [5 7 7 6 8 7 20 0]'*1000
nt=size(cashFlows,1)
nb=size(cashFlows,2)
Rates = [0.01 0.015 0.017 0.019 0.02 0.025 0.027 0.029]'

Disc= [0.99009901 0.970661749 0.950686094 0.927477246 0.90573081 0.862296866 0.829863942 0.795567442]

#Number of bonds available
nBonds = [10 100 20 30 5]'

# Declaring variables
@variable(myModel, 0<= x[b = 1:length(nBonds)] <= nBonds[b])

# Setting the objective
# The objective is to determine the cheapest set of bonds to cover liability cash flows subject to the constraint below.

@objective(myModel,Min,prices*x) 

# Add constraint
# The constraint is determined by cash flows shortfall calculated as net cash flows (bonds cash flows - liability cash flows) 
# compounded at forward rate (B), and discounted to time zero. The maximum shortfall should be less than or equal to 0.05*53.844.
# The max is linearized by setting each element in the matrix to be less than or equal to 0.05*53.844
B=[1.01 1.020024752 1.02101183 1.02502363 1.024009823 1.050370059 1.039082218 1.043109481]'
A=-(cashFlows*x-Liab_CFs)
d=cumprod(B[1:end-1])
d=[1 d[:]]'
M=tril(d./transpose(d))
@constraint(myModel, constraint1, ((M*A).*Disc)<=0.05*53.844)
# M*A is doing the following operation
# M*A=[A[1];A[1]*B[1]+A[2];(A[1]*B[1]+A[2])*B[2]+A[3];(A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4];...
#    ((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5]
#    (((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6];...
#   ((((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6])*B[6]+A[7];...
#   (((((A[1]*B[1]+A[2])*B[2]+A[3])*B[3]+A[4])*B[4]+A[5])*B[5]+A[6])*B[6]+A[7])*B[7]+A[8]]

# Printing the prepared optimization model
print(myModel)

# Solving the optimization problem
optimize!(myModel)

# Printing the optimal solutions obtained
println("Optimal Solutions:")
println("x = ", JuMP.value(x))
4

Your issue is that prices * x is an array (of length 1), and JuMP expects a scalar.

You need to use one of the following

julia> using JuMP

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
@vaModel mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, x[1:2])
2-element Array{VariableRef,1}:
 x[1]
 x[2]

julia> p = [1 2]
1×2 Array{Int64,2}:
 1  2

julia> p * x
1-element Array{GenericAffExpr{Float64,VariableRef},1}:
 x[1] + 2 x[2]

julia> @objective(model, Min, sum(p[i] * x[i] for i = 1:2))
x[1] + 2 x[2]

julia> q = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> @objective(model, Min, q' * x)
x[1] + 2 x[2]

julia> using LinearAlgebra

julia> @objective(model, Min, dot(p, x))
x[1] + 2 x[2]
5

Thank you. I have managed to fix my code but there is another point I would like to check. I have defined the constraint using a FOR loop, but I am trying to avoid the loop
and use vectorized code.

I have given my code that runs without problems at the bottom. If I replace the constraint with the following vectorized code, it results in error message. The error message suggests that it is detecting it as a quadratic expression, but I cannot see any reason why that is the case. One possible reason is that the upper triangle elements in the array M are appearing as dots or question marks rather than zero.

d=cumprod(B[1:end-1])
d=[1; d]
M1=ones(8,8).*d
M=LowerTriangular(M1./M1')
@constraint(myModel, constraint1, ((M*A).*Disc).<=0.05*53.844)

MethodError: Cannot convert an object of type GenericQuadExpr{Float64,VariableRef} to an object of type GenericAffExpr{Float64,VariableRef}
Closest candidates are:
convert(::Type{GenericAffExpr{T,V}}, !Matched::Real) where {T, V} at C:\Users\a.julia\packages\JuMP\MsUSY\src\aff_expr.jl:295
convert(::Type{GenericAffExpr{T,V}}, !Matched::V) where {T, V} at C:\Users\a.julia\packages\JuMP\MsUSY\src\aff_expr.jl:294
convert(::Type{T}, !Matched::GenericAffExpr{T,VarType} where VarType) where T at C:\Users\a.julia\packages\JuMP\MsUSY\src\aff_expr.jl:298

Stacktrace:
[1] setindex! at .\array.jl:784 [inlined]
[2] lmul!(::LowerTriangular{GenericAffExpr{Float64,VariableRef},Array{GenericAffExpr{Float64,VariableRef},2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\triangular.jl:833
[3] *(::LowerTriangular{Float64,Array{Float64,2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\triangular.jl:1874
[4] top-level scope at C:\Users\a.julia\packages\JuMP\MsUSY\src\macros.jl:390
[5] top-level scope at In[5]:38

The code below runs without errors.

using JuMP
using Cbc
using LinearAlgebra
myModel = Model(with_optimizer(Cbc.Optimizer))
prices = [99.74, 91.22, 98.71, 103.75, 97.15]
cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 102.5 5 4;4 5 0 105 104;4 105 0 0 0; 104 0 0 0 0]
Liab_CFs = [5, 7, 7, 6, 8, 7, 20, 0]*1000
nt=size(cashFlows,1)
nb=size(cashFlows,2)
Rates = [0.01, 0.015, 0.017, 0.019, 0.02, 0.025, 0.027, 0.029]
Disc= [0.99009901, 0.970661749, 0.950686094, 0.927477246, 0.90573081, 0.862296866, 0.829863942, 0.795567442]'
nBonds = [10, 100, 20, 30, 5]*1000
@variable(myModel, 0<= x[b = 1:length(nBonds)] <= nBonds[b])
@objective(myModel,Min,prices' *x) 
B=[1.01, 1.020024752, 1.02101183, 1.02502363, 1.024009823, 1.050370059, 1.039082218, 1.043109481]
A=Liab_CFs-cashFlows*x
M=zeros(length(A))
M=A[1]
for k=2:length(A)
    M(k)=M[k-1]*B[k-1]+A[k]
end
@constraint(myModel, constraint1,M.*Disc .<=0.05*53.844)
print(myModel)
optimize!(myModel)
println("Optimal Solutions:")
println("x = ", JuMP.value.(x))
6

Okay, there’s a few things going on here.

  1. The LowerTriangular thing is a bug, or at least, something isn’t getting promoted correctly in the right place. (Edit: I’ve opened an issue Issues · jump-dev/JuMP.jl · GitHub)

The current work-around is to use M = collect(LowerTriangular(...)).

  1. On for-loops vs vectorized code:

Unlike MATLAB, there is nothing wrong with using for loops in Julia. If the for loop works, keep it.

  1. On your code running without error:

M=A[1]
for k=2:length(A)
M(k)=M[k-1]*B[k-1]+A[k]
end
@constraint(myModel, constraint1,M.Disc .<=0.0553.844)

This isn’t doing what you think it’s doing. M(k) in the loop is defining a function.

7

Thanks. To fix the for loop, I changed the code to the following, but that results in “InexactError: convert…” error. Is this an issue with Julia/JuMP or I am doing something incorrectly?

M=zeros(8)
M[1] = A[1]
for k in 2:length(A)
M[k]=M[k-1]*B[k-1]+A[k]
end

InexactError: convert(Float64, -4 x[1] - 5 x[2] - 2.5 x[3] - 5 x[4] - 4 x[5] + 5000)

Stacktrace:
[1] convert at C:\Users\a.julia\packages\JuMP\MsUSY\src\aff_expr.jl:299 [inlined]
[2] setindex!(::Array{Float64,1}, ::GenericAffExpr{Float64,VariableRef}, ::Int64) at .\array.jl:782
[3] top-level scope at In[20]:18

8

The issue is that M=zeros(8) makes a Vector{Float64}, which you then try to store an affine expression into.

Use M = [AffExpr(0.0) for _ = 1:8] instead.

9

Thank you