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

Specific Domains Optimization (Mathematical)
1

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

#auxiliar function for better reading
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
2

It looks like this is an integer programming problem?

3

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

4

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

5

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

6

This is probably more in the domain of JuMP

2 Likes
7

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

8

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

image

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

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.

10

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.

11

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
12

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]
Cost(L,i) = Eval_BasicLoadCost(L, L.s[i])

#this is the struct i'm working on
 struct BasicLoad
        # (_release time_)
        r::Time  
        #(_expected ON time_)      
        e::Time       
        # (_dead line_)  
        d::Time          
        #(_start time_)
        s::StepRange{Time, Minute} 
        # discrete time step
        Δt::Time 
        #Load duration 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   
        #Load Average Power           
        P_av::Number  
        #Load Peak Power 
        P_pk::Number     

        BasicLoad() = new()

        function BasicLoad(r,e,d,Δt,C,P; allow_rand_peak = true)
		    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
13

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.

Do instead:

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

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]
5. <mark>**JuMP_VectorCost**(::Vector{Main.var"RAW_IDEAS.jl".BasicLoad})</mark>@*[Other: 13](http://localhost:1234/edit?id=7f4cfdf0-cbfd-11ed-3860-9df06b1a7405#)*
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!!!

15

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

1 Like
16

@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
17

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
18

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