# Quadratic Expression in JuMP
**URL:** https://discourse.julialang.org/t/quadratic-expression-in-jump/98601
**Category:** Optimization (Mathematical)
**Tags:** jump
**Created:** [May 10, 2023, 9:04am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601 "2023-05-10T09:04:35Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 10, 2023, 9:04am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/1 "2023-05-10T09:04:35Z")
```julia
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 10, 2023, 9:39am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/2 "2023-05-10T09:39:57Z")
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).
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 10, 2023, 9:50am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/3 "2023-05-10T09:50:55Z")
```julia
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)
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 10, 2023, 10:14am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/4 "2023-05-10T10:14:33Z")
See [Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757)
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.
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 10, 2023, 11:41am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/5 "2023-05-10T11:41:03Z")
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 11, 2023, 1:18am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/6 "2023-05-11T01:18:08Z")
> 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…
```julia
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
```Julia
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:
```julia
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
```Julia
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 11, 2023, 3:38am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/7 "2023-05-11T03:38:33Z")
You can also simplify a lot o how you’re writing the constraints. This might help:
```plaintext
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
```
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 11, 2023, 4:46am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/8 "2023-05-11T04:46:30Z")
Thank you for the simplified code. Removing instances of PgenT worked.
```julia
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
```julia
transpose (Pgen)*cov*Pgen
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 11, 2023, 5:33am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/9 "2023-05-11T05:33:37Z")
> `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
```julia
@variable(model, Pgen[bus_key, scenarios] >= 0, container = Array)
sum(Pgen[:, j]' * cov * Pgen[:, j] for j in scenarios)
```
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 11, 2023, 10:09am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/10 "2023-05-11T10:09:08Z")
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.
```julia
sum( Pgen[i, j] * cov[i, k] * Pgen[k, j]
for i in bus_key for j in scenarios for k in bus_key)
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 11, 2023, 10:29pm UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/11 "2023-05-11T22:29:00Z")
> [@Gagan\_Meena](#):
>
> a product of three matrices.
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)` ?
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 12, 2023, 4:53am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/12 "2023-05-12T04:53:41Z")
Sorry for the incomplete information. You are absolutely correct. The objective function if the sum of all the entries of the resultant matrix.
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 12, 2023, 4:59am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/13 "2023-05-12T04:59:29Z")
> [@odow](#):
>
> `sum(Array(Pgen[:, j])' * cov * Array(Pgen[:, j]) for j in scenarios)`
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
> [@odow](#):
>
> `sum(Array(Pgen[:, j])' * cov * Array(Pgen[:, k]) for j in scenarios for k in scenarios)`
Thank you very much for taking out time for this problem. Your help is highly appreciated.
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 12, 2023, 7:43am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/14 "2023-05-12T07:43:27Z")
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?
---
### Author:  [@Gagan\_Meena](https://discourse.julialang.org/u/Gagan_Meena)
#### Post date: [May 12, 2023, 7:53am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/15 "2023-05-12T07:53:57Z")
This happens only for a particular value of P=10^15, beyond which the model becomes infeasible.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [May 13, 2023, 2:28am UTC](https://discourse.julialang.org/t/quadratic-expression-in-jump/98601/16 "2023-05-13T02:28:24Z")
> [@Gagan\_Meena](#):
>
> Both of them are resulting in different values of the objective function
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
- [https://www.gurobi.com/documentation/9.5/refman/guidelines\_for\_numerical\_i.html](https://www.gurobi.com/documentation/9.5/refman/guidelines_for_numerical_i.html)
- [https://www.gurobi.com/documentation/9.5/refman/grb\_tolerances\_and\_the\_lim.html](https://www.gurobi.com/documentation/9.5/refman/grb_tolerances_and_the_lim.html)