Please help me with my optimal power flow problem

General Usage
,
1

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
α=[30.0053 30.00889 30.00741]
β=[311.669 310.333 310.833]
γ=[3213.1 3200 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(theta1im)
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)

2

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

3

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

1 Like