 # Linear programming with Jump

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)

# 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)

# 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, ((M
A).Disc)<=0.0553.844)

print(myModel)

# Solving the optimization problem

optimize!(myModel)

# Printing the optimal solutions obtained

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

The likely issue is that you need:

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

Here is the relevant documentation: https://www.juliaopt.org/JuMP.jl/stable/variables/

1 Like

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 + 91.22 x + 98.71 x + 103.75 x + 97.15 x]` is not supported by JuMP.

Stacktrace:
 error(::String) at .\error.jl:33
 set_objective(::Model, ::MathOptInterface.OptimizationSense, ::Array{GenericAffExpr{Float64,VariableRef},1}) at C:\Users\a.julia\packages\JuMP\MsUSY\src\objective.jl:101
 top-level scope at C:\Users\a.julia\packages\JuMP\MsUSY\src\macros.jl:979
 top-level scope at In: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)

# 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;A*B+A;(A*B+A)*B+A;(A*B+A)*B+A)*B+A;...
#    ((A*B+A)*B+A)*B+A)*B+A
#    (((A*B+A)*B+A)*B+A)*B+A)*B+A;...
#   ((((A*B+A)*B+A)*B+A)*B+A)*B+A)*B+A;...
#   (((((A*B+A)*B+A)*B+A)*B+A)*B+A)*B+A)*B+A]

# 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))
``````

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
x

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

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

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

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

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

julia> using LinearAlgebra

julia> @objective(model, Min, dot(p, x))
x + 2 x
``````

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:
 setindex! at .\array.jl:784 [inlined]
 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
 *(::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
 top-level scope at C:\Users\a.julia\packages\JuMP\MsUSY\src\macros.jl:390
 top-level scope at In: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
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))
``````

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 https://github.com/JuliaOpt/JuMP.jl/issues/2157.)

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
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.

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 = A
for k in 2:length(A)
M[k]=M[k-1]*B[k-1]+A[k]
end

InexactError: convert(Float64, -4 x - 5 x - 2.5 x - 5 x - 4 x + 5000)

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

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.

Thank you