Quadratic Expression in JuMP

Specific Domains Optimization (Mathematical)
1 
  ESS_cost = CRF*(CapES*max_Eb*100 + CapPS*max_Pb*100) + FOM*max_Pb*100 + VOM*sum(PD[s]*dur[s] for s in scenar)*100*365
    cost_genvar = sum(PgenT[j,i]*cov[i,k]*Pgen[k,j] for i in bus_key for j in scenar for k in bus_key)
    
    if p == 0
        @objective(ESS, Min, ESS_cost)
    else
        @objective(ESS, Min, ESS_cost + p*cost_genvar)
    end

Is this the right way to write a quadratic expression for a QP problem? My problem formation shows that Q is not PSD. Is that because of the incorrect implementation of quadratic expression? This cannot be because of data as I have used the same data for writing the same code in python and it worked.

2

Hi there,

Please provide a fully reproducible example. It’s difficult to say what’s going on without one. What are the decisions variables and what are constants?

If you’re expecting your problem to be convex, then there might be a mistake. But expressions like x * y are not convex.

If you know your problem is not convex, some solvers support non-convex quadratics. You could try Gurobi (if there are integer variables) or Ipopt (if everything is continuous).

3 
    ESS = Model(Gurobi.Optimizer)

    # Variables
    Pgen = @variable(ESS, [bus_key, scenar], lower_bound=0, base_name ="Pgen")
    PD = @variable(ESS, [scenar], lower_bound=0, base_name="PD")
    PC = @variable(ESS, [scenar], lower_bound=0, base_name="PC")
    Eb = @variable(ESS, [scenar], lower_bound=0, base_name="Eb")
    max_Pb = @variable(ESS, lower_bound=0, base_name="max_Pb")
    max_Eb = @variable(ESS, lower_bound=0, base_name="max_Eb")
    PgenT = @variable(ESS, [scenar, bus_key], lower_bound=0, base_name ="PgenT")
    #cost_genvar = @variable(ESS)

    ## constraints
    Pgenlim = @constraint(ESS, [n=bus_key,s=sc], Pgen[n,s]<=PG[n,s], base_name="Pgenlim")
    kcl = @constraint(ESS, [s=sc], sum(Pgen[n,s] for n in bus_key) + PD[s] - PC[s] - PL[s] == 0, base_name="kcl")
    
    # Limits on Charging and Discharging Power, Energy, and Pmax, Emax
    PDmax_up_lim =@constraint(ESS, [s=sc], PD[s] <= ESS_stat*max_Pb, base_name="PDmax_up_lim")
    PCmax_up_lim =@constraint(ESS, [s=sc], PC[s] <= ESS_stat*max_Pb, base_name="PCmax_up_lim")
    
    # Limits on Emax
    Eb_up_lim = @constraint(ESS, [s=sc], Eb[s] <= ESS_stat*0.8*max_Eb, base_name="Eb_up_lim")
    Eb_low_lim = @constraint(ESS, [s=sc], Eb[s] >= ESS_stat*0.2*max_Eb, base_name="Eb_low_lim")
    Pbmax_up_lim = @constraint(ESS, max_Pb <= 10^6, base_name="Pbmax_up_lim")

    # Limit on Discharge Duration
    Ebmax_up_lim = @constraint(ESS, max_Eb <= DSmax*max_Pb/eta, base_name="Ebmax_up_lim")
    Ebmax_low_lim = @constraint(ESS, max_Eb >= DSmin*max_Pb/eta, base_name="Ebmax_low_lim")
    
    # Energy Balance constraint

    Energy_balanceall = @constraint(ESS,  [s=scenar_except_end], Eb[s+1] - Eb[s] + (PD[s]/eta - PC[s]*eta)*dur[s] == 0, base_name="Energy_balanceall") 
    Energy_balanceend = @constraint(ESS,  [n_scen], Eb[1] - Eb[n_scen] + (PD[n_scen]/eta - PC[n_scen]*eta)*dur[n_scen] == 0, base_name="Energy_balanceend")

    println(Energy_balanceend)
    ## 
    CRF=0.1130
    for i in 1:n_bus
        for j in 1:n_scen
            @constraint(ESS, Pgen[i,j]==PgenT[j,i])
        end
    end

    # cost_ESS
    ESS_cost = CRF*(CapES*max_Eb*100 + CapPS*max_Pb*100) + FOM*max_Pb*100 + VOM*sum(PD[s]*dur[s] for s in scenar)*100*365
    cost_genvar = sum(PgenT[j,i]*cov[i,k]*Pgen[k,j] for i in bus_key for j in scenar for k in bus_key)
    
    if p == 0
        @objective(ESS, Min, ESS_cost)
    else
        @objective(ESS, Min, ESS_cost + p*cost_genvar)
    end
    # solve
    set_optimizer_attribute(ESS, "OutputFlag", 1)
    set_optimizer_attribute(model, "FeasibilityTol", 1e-8)
    println(cost_genvar)
4

See Please read: make it easier to help you

I cant copy and paste your code so i cant run it to find the problem.

Whats the full text of the error message?

Try simplifying the problem. Does it solve without the objective? Do you get the same error if you delete all constraints?

Try to narrow down whats causing the issue and then post an example that someone can copy paste and trigger the same issue.

5

Sorry for the inconvenience. However, I am able to copy the code.

This is the full error message.

Gurobi Error 10020: Objective Q not PSD (diagonal adjustment of 1.4e+07 would be required). Set NonConvex parameter to 2 to solve model.
Stacktrace:
[1] _check_ret
@ C:\Users\Admin.julia\packages\Gurobi\EKa6j\src\MOI_wrapper\MOI_wrapper.jl:375 [inlined]
[2] _check_ret_GRBoptimize(model::Gurobi.Optimizer)
@ Gurobi C:\Users\Admin.julia\packages\Gurobi\EKa6j\src\MOI_wrapper\MOI_wrapper.jl:392
[3] optimize!(model::Gurobi.Optimizer)
@ Gurobi C:\Users\Admin.julia\packages\Gurobi\EKa6j\src\MOI_wrapper\MOI_wrapper.jl:2674
[4] optimize!
@ C:\Users\Admin.julia\packages\MathOptInterface\8f6oN\src\Bridges\bridge_optimizer.jl:376 [inlined]
[5] optimize!
@ C:\Users\Admin.julia\packages\MathOptInterface\8f6oN\src\MathOptInterface.jl:85 [inlined]
[6] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
@ MathOptInterface.Utilities C:\Users\Admin.julia\packages\MathOptInterface\8f6oN\src\Utilities\cachingoptimizer.jl:316
[7] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ JuMP C:\Users\Admin.julia\packages\JuMP\AKvOr\src\optimizer_interface.jl:440
[8] optimize!
@ C:\Users\Admin.julia\packages\JuMP\AKvOr\src\optimizer_interface.jl:410 [inlined]
[9] ESS_Sizing(L::Ld, E::St, sc::Vector{Int64}, cov::Matrix{Float64}, p::Float64)
@ Main d:\Windows PC\IITB\Sem5\Optimization with Julia\ESS_Sizing.jl:227
[10] top-level scope
@ .\timing.jl:273

Yes, the problem is solvable without the objective function. It doesn’t result in any error.

The second component of the objective function which is cost_genvar is causing the error.

6

I cant copy and paste your code so i cant run it to find the problem.
However, I am able to copy the code.

Sure, I can copy the code. But you haven’t defined any of the data so it won’t actually run.

But oh, I see what you’ve done…

PgenT = @variable(ESS, [scenar, bus_key], lower_bound=0, base_name ="PgenT")

Is this meant to be the transpose of Pgen???

I guess so, because you have

    for i in 1:n_bus
        for j in 1:n_scen
            @constraint(ESS, Pgen[i,j]==PgenT[j,i])
        end
    end

Just delete all occurances of PgenT and use:

cost_genvar = sum(Pgen[i,j]*cov[i,k]*Pgen[k,j] for i in bus_key for j in scenar for k in bus_key)

If you have something like

model = Model()
@variable(model, x)
@variable(model, y)
@constraint(model, x == y)
@objective(model, Min, x * y)

Gurobi will complain that the problem is non-convex, even though you as the author know that it is because of the x == y constraint. The objective needs to be convex excluding any constraints.

2 Likes
7

You can also simplify a lot o how you’re writing the constraints. This might help:

using JuMP, Gurobi
CRF = 0.1130
model = Model(Gurobi.Optimizer)
set_attribute(model, "OutputFlag", 1)
set_attribute(model, "FeasibilityTol", 1e-8)
@variables(model, begin
    Pgen[bus_key, scenarios] >= 0
    PD[scenarios] >= 0
    PC[scenarios] >= 0
    Eb[scenarios] >= 0
    max_Pb >= 0
    0 <= max_Eb <= 1e6
end)
@constraints(model, begin
    [n=bus_key, s=scenarios], Pgen[n, s] <= PG[n, s]
    [s=scenarios], sum(Pgen[:, s]) + PD[s] - PC[s] - PL[s] == 0
    [s=scenarios], PD[s] <= model_stat * max_Pb
    [s=scenarios], PC[s] <= model_stat * max_Pb
    [s=scenarios], Eb[s] <= 0.8 * model_stat * max_Eb
    [s=scenarios], 0.2 * model_stat * max_Eb <= Eb[s]
    max_Eb <= DSmax * max_Pb / eta
    DSmin * max_Pb / eta <= max_Eb
    [s=scenarios[1:end-1]], Eb[s+1] - Eb[s] + (PD[s]/eta - PC[s]*eta)*dur[s] == 0
    [s=scenarios[end]], Eb[1] - Eb[s] + (PD[s]/eta - PC[s]*eta)*dur[s] == 0
end)
@expressions(model, begin
    model_cost, CRF * 100 * (CapES * max_Eb + CapPS * max_Pb) + 
                FOM * 100 * max_Pb + 
                VOM * 100 * 365 * sum(PD[s] * dur[s] for s in scenarios)
    cost_genvar, sum(
                    Pgen[i, j] * cov[i, k] * Pgen[k, j]
                    for i in bus_key for j in scenarios for k in bus_key
                 )
end)
if p == 0
    @objective(model, Min, model_cost)
else
    @objective(model, Min, model_cost + p * cost_genvar)
end
1 Like
8

Thank you for the simplified code. Removing instances of PgenT worked.

cost_genvar, sum(
                    Pgen[i, j] * cov[i, k] * Pgen[k, j]
                    for i in bus_key for j in scenarios for k in bus_key
                 )

Here, the above product will simply square all the terms of the matrix Pgen. How to achieve the following

transpose (Pgen)*cov*Pgen
9

transpose (Pgen)*cov*Pgen

Pgen is a matrix. I don’t know if it makes sense to write that.

Did you meant instead sum(Array(Pgen[:, j])' * cov * Array(Pgen[:, j]) for j in scenarios)?

Alternatively, if bus_key and scenarios are both 1-based integer indices like 1, 2, 3, ..., N, then you could use

@variable(model, Pgen[bus_key, scenarios] >= 0, container = Array)
sum(Pgen[:, j]' * cov * Pgen[:, j] for j in scenarios)
10

Since Pgen is a matrix of size [bus_key, scenarios]
cov is a matrix of size [bus_key, bus_key]
I want to write the quadratic expression as a product of three matrices.
1). transpose of Pgen
2). cov
3). Pgen

cov is a constant matrix

using the product mentioned below I am not getting all the terms in the product of these three matrices.

sum( Pgen[i, j] * cov[i, k] * Pgen[k, j]
     for i in bus_key for j in scenarios for k in bus_key)
11

The result is a matrix of size [scenarios, scenarios]. That cannot be your objective function, so I think you want something else. My guess is:

Did you meant instead sum(Array(Pgen[:, j])' * cov * Array(Pgen[:, j]) for j in scenarios) ?

12

Sorry for the incomplete information. You are absolutely correct. The objective function if the sum of all the entries of the resultant matrix.

13

But it is not just j the jth column will be multiplied with all the columns. Thanks a lot now I got this. It should be

Thank you very much for taking out time for this problem. Your help is highly appreciated.

2 Likes
14

I used the compact formulation that you mentioned by suing @variables and @constraints.
I solved the problem using the formulation that you gave and the formulation that I used. Both of them are resulting in different values of the objective function. This should not happen.
What might be the possible cause?

15

This happens only for a particular value of P=10^15, beyond which the model becomes infeasible.

16

Then they’re not equivalent. but it’s hard to know the exact difference because you haven’t provided a reproducible example. At a guess, I think your for j in scenarios for k in scenarios is wrong. Did you try my original suggestion?

This happens only for a particular value of P=10^15, beyond which the model becomes infeasible.

Using such a large value of P will run into numerical error and you should not expect the solver to find a solution (it probably printed a bunch of warnings?).

See