# Please help me with my optimal power flow problem

**URL:** <https://discourse.julialang.org/t/please-help-me-with-my-optimal-power-flow-problem/31211>\
**Category:** General Usage\
**Tags:** question, homework\
**Created:** [November 18, 2019, 7:16am UTC](https://discourse.julialang.org/t/please-help-me-with-my-optimal-power-flow-problem/31211 "2019-11-18T07:16:09Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![varun\_yaramasu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/varun_yaramasu/32/10336_2.png) [@varun\_yaramasu](https://discourse.julialang.org/u/varun_yaramasu)\
**Post date:** [November 18, 2019, 7:16am UTC](https://discourse.julialang.org/t/please-help-me-with-my-optimal-power-flow-problem/31211/1 "2019-11-18T07:16:09Z")

</div>

using JuMP,Ipopt

m=Model(with\_optimizer(Ipopt.Optimizer))

# from to r x b Max\_MVA(rate)

branch=[ 1 2 0.1 0.2 0.04 60  
1 4 0.05 0.2 0.04 90  
1 5 0.08 0.3 0.06 50  
2 3 0.05 0.25 0.06 40  
2 4 0.05 0.1 0.02 70  
2 5 0.1 0.3 0.04 50  
2 6 0.07 0.2 0.05 50  
3 5 0.12 0.26 0.05 50  
3 6 0.02 0.1 0.02 80  
4 5 0.2 0.4 0.08 30  
5 6 0.1 0.3 0.06 30  
];

# bus bus type Pd Qd Pg Qg theta

bus=[ 1 3 0 0 110 0 0  
2 2 0 0 80 0 0  
3 2 0 0 50 0 0  
4 1 50 10 0 0 0  
5 1 100 15 0 0 0  
6 1 100 15 0 0 0  
];

n\_b=length(bus[:,1])  
ng=3;#generation buses+slack bus  
nl=3;#load buses

Ybus =[ 4.0063-11.7479im -2.0000+4.0000im 0.0000+0.0000im -1.1765+4.7059im -0.8299+3.1120im 0.0000+0.0000im  
-2.0000+4.0000im 9.3283-23.1955im -0.7692+3.8462im -4.0000+8.0000im -1.0000+3.0000im -1.5590+4.4543im  
0.0000+0.0000im -0.7692+3.8462im 4.1557-16.5673im 0.0000+0.0000im -1.4634+3.1707im -1.9231+9.6154im  
-1.1765+4.7059im -4.0000+8.0000im 0.0000+0.0000im 6.1765-14.6359im -1.0000+2.0000im 0.0000+0.0000im  
-0.8299+3.1120im -1.0000+3.0000im -1.4634+3.1707im -1.0000+2.0000im 5.2933-14.1378im -1.0000+3.0000im  
0.0000+0.0000im -1.5590+4.4543im -1.9231+9.6154im 0.0000+0.0000im -1.0000+3.0000im 4.4821-17.0047im  
]

@variable(m,PG[1:ng])  
@variable(m,QG[1:ng])  
@variable(m,theta[1:(n\_b-1)])  
@variable(m,V[1:n\_b])

#D=250  
α=[3_0.0053 3_0.00889 3_0.00741]  
β=[3_11.669 3_10.333 3_10.833]  
γ=[3_213.1 3_200 3\*240]

P\_min=[50 37.5 45]’  
P\_max=[200 150 180]’

Q\_min=[-100 -100 -100]’  
Q\_max=[150 150 120]’

V\_min=[0.9 0.9 0.9 0.9 0.9 0.9]’  
V\_max=[1.1 1.1 1.1 1.1 1.1 1.1]’

theta\_min=[-20 -20 -20 -20 -20]‘\*pi/180  
theta\_max=[20 20 20 20 20]’\*pi/180

LB=[theta\_min; V\_min; P\_min; Q\_min]  
UB=[theta\_max; V\_max; P\_max; Q\_max]

Pg=[PG[1] PG[2] PG[3] 0 0 0]‘/100  
Qg=[QG[1] QG[2] QG[3] 0 0 0]’/100  
Pd=[0 0 0 50 100 100]‘/100 #%Active Power demand at each bus  
Qd=[0 0 0 10 15 15]’/100 #%Reactive Power demand at each bus

P=Pg-Pd #Active power calculated at ith bus  
Q=Qg-Qd # Reactive power calculated ith bus

@variable(m,Paux)  
@constraint(m,Paux .==P)

@variable(m,Qaux)  
@constraint(m,Qaux .==Q)

theta=[0;theta]

function V\_bus(V,theta)  
V._exp(theta_1im)  
end

function S\_bus(V\_bus,Ybus)  
V\_bus(V,theta).\*(conj(Ybus.\*V\_bus(V,theta)))  
end

function P\_bus(S\_bus)  
real(S\_bus(V\_bus(V,theta),Ybus))  
end

function Q\_bus(S\_bus)  
imag(S\_bus(V\_bus(V,theta),Ybus))  
end

@objective(m,Min,sum(α[i]\*PG[i]^2+β[i]\*PG[i]+γ[i] for i in 1:3))

@constraint(m,Paux-P\_bus(S\_bus(V\_bus(V,theta),Ybus)) ==0)  
@constraint(m,Qaux-Q\_bus(S\_bus(V\_bus(V,theta),Ybus)) ==0)

@constraint(m,theta[1:5] .\<=UB[1:5])  
@constraint(m,theta[1:5] .\>=LB[1:5])

@constraint(m,V[1:n\_b] .\<=UB[6:11])  
@constraint(m,V[1:n\_b] .\>=LB[6:11])

@constraint(m,PG[1:ng] .\<=UB[12:14])  
@constraint(m,PG[1:ng] .\>=LB[12:14])

@constraint(m,QG[1:ng] .\<=UB[15:17])  
@constraint(m,QG[1:ng] .\>=LB[15:17])

println(“The optimization problem to be solved is:”)  
print(m)  
JuMP.optimize!(m)  
optimvalue=JuMP.objective\_value(m)  
PG=JuMP.value.(PG)  
QG=JuMP.value.(QG)  
V=JuMP.value.(V)  
theta=JuMP.value.(theta)

println(“Objective value:”,optimvalue)  
println(“PG1=:”,PG[1])  
println(“PG2=:”,PG[2])  
println(“PG3=:”,PG[3])

println(“QG1=:”,QG[1])  
println(“QG2=:”,QG[2])  
println(“QG3=:”,QG[3])

println(“V1=:”,V[1])  
println(“V2=:”,V[2])  
println(“V3=:”,V[3])  
println(“V4=:”,V[4])  
println(“V5=:”,V[5])  
println(“V6=:”,V[6])

println(“theta2=:”,theta[1])  
println(“theta3=:”,theta[2])  
println(“theta4=:”,theta[3])  
println(“theta5=:”,theta[4])  
println(“theta6=:”,theta[5])

#PG3=JuMP.value(PG[i] for i in 3)  
#PG2=JuMP.value(PG[i])  
#PG3=JuMP.value(PG3)

#println(“PG2=:”,PG2)  
#println(“PG3=:”,PG3)

---

<div class="post-metadata">

**Author:** ![varun\_yaramasu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/varun_yaramasu/32/10336_2.png) [@varun\_yaramasu](https://discourse.julialang.org/u/varun_yaramasu)\
**Post date:** [November 18, 2019, 7:17am UTC](https://discourse.julialang.org/t/please-help-me-with-my-optimal-power-flow-problem/31211/2 "2019-11-18T07:17:14Z")

</div>

julia\>  
ERROR: LoadError: InexactError: Float64(im)  
Stacktrace:  
[1] Type at .\complex.jl:37 [inlined]  
[2] convert at .\number.jl:7 [inlined]  
[3] \_float at C:\Users\user.juliapro\JuliaPro\_v1.2.0-1\packages\JuMP\iGamg\src\operators.jl:13 [inlined]  
[4] \* at C:\Users\user.juliapro\JuliaPro\_v1.2.0-1\packages\JuMP\iGamg\src\operators.jl:42 [inlined]  
[5] \* at C:\Users\user.juliapro\JuliaPro\_v1.2.0-1\packages\JuMP\iGamg\src\operators.jl:123 [inlined]  
[6] \_broadcast\_getindex\_evalf at .\broadcast.jl:625 [inlined]  
[7] \_broadcast\_getindex at .\broadcast.jl:598 [inlined]  
[8] getindex at .\broadcast.jl:558 [inlined]  
[10] macro expansion at .\simdloop.jl:77 [inlined]  
[11] copyto! at .\broadcast.jl:887 [inlined]  
[12] copyto! at .\broadcast.jl:842 [inlined]  
[13] copy at .\broadcast.jl:818 [inlined]  
[14] materialize at .\broadcast.jl:798 [inlined]  
[15] broadcast(::typeof(\*), ::Array{GenericAffExpr{Float64,VariableRef},1}, ::Complex{Bool}) at .\broadcast.jl:752  
[16] \* at .\arraymath.jl:55 [inlined]  
[17] V\_bus(::Array{VariableRef,1}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at C:\Users\user\OneDrive - The University of Colorado Denver\Power System Operation and control\Homework3\question1.jl:98  
[18] top-level scope at C:\Users\user.juliapro\JuliaPro\_v1.2.0-1\packages\JuMP\iGamg\src\macros.jl:390  
in expression starting at C:\Users\user\OneDrive - The University of Colorado Denver\Power System Operation and control\Homework3\question1.jl:144

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [November 18, 2019, 9:02am UTC](https://discourse.julialang.org/t/please-help-me-with-my-optimal-power-flow-problem/31211/3 "2019-11-18T09:02:25Z")

</div>

Welcome to the Julia Discourse! You’ve posted quite a bit of code for people to go through, you might want to consider checking out [Please read: make it easier to help you](https://discourse.julialang.org/t/psa-make-it-easier-to-help-you/14757) and follow the tips there to amend your question to increase the chance of getting help quickly.

With that said, your code errors because you are trying to convert a complex number into a floating point number somewhere, which doesn’t work. An MWE for this would be:

```julia
julia> a = 1.0 + 1im
1.0 + 1.0im

julia> typeof(a)
Complex{Float64}

julia> Float64(a)
ERROR: InexactError: Float64(1.0 + 1.0im)
Stacktrace:
 [1] Float64(::Complex{Float64}) at .\complex.jl:37
 [2] top-level scope at REPL[3]:100  

```

You would therefore have to think about where this conversion happens and why.

Also, while it is fine to ask for help with your homework here on Discourse, you might want to have a look at our informal homework policy here: [Homework policy - #17 by Tamas\_Papp](https://discourse.julialang.org/t/homework-policy/2342/17)