# Help on model a problem/function to use with Optimization.jl

Hello to all.

I’m trying to use Optimization.jl to solve a problem, but a I dont know how to define the problem in a way the package needs it. In my problem I have an Array of a Struct with several parameters and for each element of the Array I need the proper index of a field that is a StepRange type.

I’ll write here a smaller version of my struct and equation. Any help or hint is really appreciated, because I dont know even how to start it. Thanks in advance!!!

``````#Struct with problem values - needed an Array of it
struct L
a::Float64
b::Float64
s::StepRange{Time,Minute}
t1::Time
t2::Time
t3::Time

L() = new()
function L(a,b,t1,t2,t3)
s = t1:t2:t3
new(a,b,s,t1,t2,t3)
end
end

#GOAL is find all values of x that minimize this function
function Eval_EQ_SUM(H::Array{L,1},x::Array{Int64,1})
if length(H) ≠ length(x)
@error "H and x should have the same length"
return nothing
end
S = Vector{Float64}(undef,length(H))
for i in eachindex(H)
S[i] = H[i].a*EQ(H[i].s[x[i]]), H[i].b)
end
return sum(S)
end

function EQ(t::time, b::Float64)
if t < Time(2)
F = 2*b
elseif t > Time(4) && t < Time(12)
F = sqrt(b)
else
F = b^3
end
return F
end
``````

It looks like this is an integer programming problem?

Yes, values of x[i] are index to a range H[i].s and its values should be restricted to bounds of the vector H[i].s.

I’m not an expert in optimization at all but i think using integer variables should be easy to solve the problem. Past colleagues used Genetic Algorithm and PSO to solve MIN the equation, but the programing interface was Matlab/Python. Thats why Im looking into Optimization.jl

To me it looks more like a constraint programming problem, cause the syntax `s[x[i]]` is not natively supported by integer programming solvers

Well, as I said, I’m not an expert in optimization.

I’m just thinking, ‘x’ is integer so it should be an integer programming problem.

If it could help, I can replace `H[i].s[x[i]]` for a variable `t[i], typeof(t) = Vector{Time}`

My original problem is a schedule type. I should position the task within an window of time (t1:t3) in a way it minimize a not continuous function based on the schedule point

This is probably more in the domain of JuMP

2 Likes

@rodolforbcoutinho could you point me to the mathematical formulation of your optimization problem?

Here are some link for the papers I’m trying to reproduce the base model.

The image above also have the main equation I’m trying to work with. As you could see its a bit more complicated system than the one I wrote in the initial post. But I’m sure if I could model the function `Eval_EQ_SUM`, I could model the complete system.

I found a way to eliminate the inner Sum, but to do that I should call a function with some if/else logic. Dont see how could i put it in a single statemente like the one in examples (2x + 3y)

I’ll do some tests now with JuMP.jl. Its a way more clear in its tutorial “who goes where”.

In/with JuMP.jl, can I do something like this ???

``````function f(a)
if a > 10
return 2*a
elseif a < 4
return 3*a
else
return a^2
end
end

@objective(model, Min, f(x) + 20y)
``````

In/with JuMP.jl, can I do something like this ???

It depends on the rest of the problem structure. But that objective function probably doesn’t make sense, because it is discontinuous.

If you mean something like your `EQ` function, is `t` a decision variable? If it is not, then yes.

I’d start by providing a reproducible example of what you’re trying to do. Use the `f(a)` syntax for now, even if JuMP complains and gives an error. Once things are started, people can probably point you in the right direction.

The function I must handle its discontinuous. As I said before its something realy like f(a) in the last post. It allow me to cut of about 80% memory allocation and its about 10x faster than the original formulation. So I’ll keep trying on.

Also `vector t` is my decision variable.

But I’ll try to do like you said, and bring on the real problem.

I think you might want to try actually modeling the variable `u` in your JuMP problem, and using an integer programming solver like HiGHS

2 Likes

Hello again!
Finally I have some real working code to share!

But First what didn’t worked at all

``````function JuMP_Cost_alfa(L::BasicLoad)

model = Model(HiGHS.Optimizer)
u = length(L.s)
@variable(model, 1 < x < u, Int)
@objective(model, Min, Cost(L,x))
optimize!(model)

return value.(x)
end
``````

Now the code that works

``````function JuMP_Cost(L::BasicLoad)

model = Model(HiGHS.Optimizer)
u = length(L.s)
@variable(model, x[1:u], Bin)
@constraint(model,sum(x)==1)
@objective(model, Min, sum(Cost(L,i)*x[i] for i in 1:u) )
optimize!(model)

A = value.(x)
ix = findfirst(x->x==1.0,A)
return ix
end
``````

Some Complementar code

``````#this is discontinuous function with lots of logic constraints but is working fine
#just to point the use of index i inside L.s[i]
``````
``````#this is the struct i'm working on
# (_release time_)
r::Time
#(_expected ON time_)
e::Time
d::Time
#(_start time_)
s::StepRange{Time, Minute}
# discrete time step
Δt::Time
C::Time
#discrete time duration
CΔ::Int64
#Power over time model
# OBS: typeof(P) == Array{Number,1} OR Function t -> f(t)
# OBS: length(P{Array}) = CΔ
P
P_av::Number
P_pk::Number

CΔ = div(C.instant,Δt.instant)

s = r:Minute(Δt):(Time(d-C))
if typeof(P) <: Function
V = P.(1:CΔ)
elseif typeof(P) <: Array
if length(P) == CΔ
V = P
else
@error "Unable to create new Load. P::Array{Number,1} must be compatible with Δt size"
return nothing
end
else
@error "Unable to create new Load. P must be a Function t -> f(t) or Array{Number,1} type"
return nothing
end

P_av = mean(V)
P_pk = maximum(V)

if (P_pk == P_av)&&(allow_rand_peak)
P_pk = P_av*rand(1.1:0.1:1.5)
end
new(r,e,d,s,Δt,C,CΔ,P,P_av,P_pk)
end
end
``````

Yes, this approach looks good. You can’t use user-defined functions like `Cost(L, x)` with HiGHS. Everything must be explicitly written as a MILP.

You should be careful with:

``````A = value.(x)
ix = findfirst(x->x==1.0,A)
``````

HiGHS uses numerical tolerances, so you’re not guaranteed to return a value `==1.0`.

``````A = value.(x)
ix = findfirst(x -> x > 0.5, A)
# or
ix = findfirst(xi -> round(Int, value(xi)) == 1, x)
``````
1 Like

Now I’m trying to solve the same problem, but for an Array of the struct BasicLoad

Here the last code, it has a bug i couldn’t find.

``````function JuMP_VectorCost(H::Vector{BasicLoad})

model = Model(HiGHS.Optimizer)
n = length(H)
u = Vector{Int64}(undef,n)
for i in 1:n
u[i] = length(H[i].s)
end

@variable(model, x[i = n, j = u[i]], Bin)
@constraint(model, [i = n, j = u[i]], sum(x[i,:])==1)

@objective(model, Min, sum(Cost(H[i],j)*x[i,j] for i in 1:n, j in 1:u[i]) )
optimize!(model)

return value.(x)
end
``````

This is the error message

``````KeyError: key (1, 1) not found

1. <mark>**getindex**</mark>@*dict.jl:498* [inlined]
2. <mark>**getindex**(::JuMP.Containers.SparseAxisArray{JuMP.VariableRef, 2, Tuple{Int64, Int64}}, ::Int64, ::Int64)</mark>@*SparseAxisArray.jl:128*
3. <mark>**macro expansion**</mark>@*[Other: 322](http://localhost:1234/edit?id=7f4cfdf0-cbfd-11ed-3860-9df06b1a7405#)* [inlined]
4. <mark>**macro expansion**</mark>@*[Other: 1297](http://localhost:1234/edit?id=7f4cfdf0-cbfd-11ed-3860-9df06b1a7405#)* [inlined]
6. <mark>**top-level scope**</mark>@*[Local: 1](http://localhost:1234/edit?id=7f4cfdf0-cbfd-11ed-3860-9df06b1a7405#)* [inlined]
``````

Any help on the `JuMP_VectorCost function` would be nice! Thanks in advance!!!

Yeah thats a better idea indeed! The other code there was a “late night” solution kind.

1 Like

`@variable(model, x[i = n, j = u[i]], Bin)`

You probably need `@variable(model, x[i = 1:n, j = 1:u[i]], Bin)`

1 Like

It works nice now!!! I could even do it multiobjective with JuMP and HiGHS.
In full problem I should minimize COST and maximize COMFORT
Bellow the final solution for the record

``````function JuMP_VectorMulti(H::Vector{Raw.BasicLoad})

model = Model()
set_optimizer(model, () -> MOA.Optimizer(HiGHS.Optimizer))

n = length(H)
u = Vector{Int64}(undef,n)
for i in 1:n
u[i] = length(H[i].s)
end

@variable(model, x[i = 1:n, j = 1:u[i]], Bin)
@constraint(model, [i = 1:n, j = 1:u[i]], sum(x[i,:])==1)

#MIN
@expression(model, Cost_expr, sum(Cost(H[i],j)*x[i,j] for i in 1:n, j in 1:u[i]))
#MAX
@expression(model, Comf_expr, sum(Comf(H[i],j)*x[i,j] for i in 1:n, j in 1:u[i]))

@objective(model, Max, [-Cost_expr, Comf_expr])
optimize!(model)
solution_summary(model)

return value.(x)
end
``````

Should I change the title or add some other mark? maybe this solution could be usefull to others.

1 Like

I think this should be `@constraint(model, [i = 1:n], sum(x[i, :]) == 1)`.

You could also try different algorithms in MOA that return a set of solutions, not just a single point. See Multi-objective knapsack · JuMP.

1 Like