Hello,
I apologize if my question is dumb in advance. Flowing is my code, when I run it, it returns
" no method matching coeftype(::Array{JuMP.GenericAffExpr{Float64,JuMP.Variable},1})" .I don’s really understand the error message. Could any one tells me what does it mean?
using JuMP
using Clp
using Ipopt
include("acopf_input.jl")
path = pwd()"/case.dat"
casefilename, refbus, loadscale, thermalLimitscale, extras, mvaBase = readconfig(path)
generators, generatorlist, buses, lines = readcase(casefilename, extras, loadscale, thermalLimitscale, mvaBase)
numbuses = length(buses)
numlines = length(lines)
numgenerators = length(generators)
vold=[1.05 0;1.06 -0.03;1.04 0.149;1.04 0.093;1.015 0.1685;1.0331 0.1534;1.025 0.149;1.045 0.073;1.009 0.156]
#i0 is branch current
iold=[1.5 -0.7;0.78 -0.38;-0.12 0.08;0.06 -0.43;-0.05 -0.35;-1.04 0.144;-1.58 0.52;0.54 -0.37;-0.71 0.33]
ibus_opt=[-1.5 0.7;-1.58 0.52;-0.06 0.43;0 0;0.9 -0.46;0 0;0.99 -0.500 ;0 0;1.25 -0.72]
pgold=zeros(Float64,numbuses)
maxitr=30
h = 1
vsqold= zeros(Float64,numbuses)
isqold= zeros(Float64,numlines)
μv=zeros(Float64,numbuses) #!
μi=zeros(Float64,numlines) #!
P_tay=zeros(Float64,numbuses)
Q_tay=zeros(Float64,numbuses)
P_bilinear = zeros(Float64, numbuses)
Q_bilinear = zeros(Float64, numbuses)
βr = zeros(Float64, numbuses)
βj = zeros(Float64, numbuses)
αr = zeros(Float64, numbuses)
αj = zeros(Float64, numbuses)
φpn = zeros(Float64, numbuses)
φqn = zeros(Float64, numbuses)
stepvr = zeros(Float64, numbuses)
stepvj = zeros(Float64, numbuses)
stepir = zeros(Float64, numlines)
stepij= zeros(Float64, numlines)
δp=ones(Float64,numbuses)
δq=ones(Float64,numbuses)
δp_tot=numbuses
δq_tot=numbuses
while((h<30)&&((maximum(δp)>1e-4)&&(maximum(δq)>1e-4))&&((δp_tot>5*1e-4)&&(δq_tot>5*1e-4)))
function solve_taylor(thermalLimitscale, generators, buses, lines)
numbuses = length(buses)
numlines = length(lines)
numgenerators = length(generators)
m=Model(solver=IpoptSolver())
@variable(m,ir[1:numbuses]) #real current injected at node
@variable(m,ij[1:numbuses]) #immaginary current injected at node
@variable(m,irij[1:numlines]) #real current flow from i to j
@variable(m,ijij[1:numlines]) #immaginary current flow from i to j
@variable(m, gp[1:numgenerators])#Active power generated at certain generator
@variable(m, gq[1:numgenerators])#Reactive power generated
@variable(m, vr[1:numbuses] >=0)#Bus voltage real part
@variable(m, vj[1:numbuses])#Bus voltage imanginary part
@variable(m,vsqnew[1:numbuses] >= 0)
@variable(m,isqnew[1:numlines] >= 0)
@variable(m,gp_tot[1:numbuses])
@variable(m,gq_tot[1:numbuses])
@objective(m, Min, sum{generators[i].pi1*gp[i]^2+generators[i].pi2*gp[i]+generators[i].pi3, i=1:numgenerators})
for j = 1:numlines
con=@constraint(m,irij[j] == lines[j].γ* vr[lines[j].head] -lines[j].β * vj[lines[j].head] - lines[j].γ* vj[lines[j].tail] + lines[j].β*vj[lines[j].tail] )
con=@constraint(m,ijij[j] == lines[j].γ* vj[lines[j].head] + (lines[j].β-lines[j].b_charge/2) * vr[lines[j].head] - lines[j].γ* vj[lines[j].tail] -lines[j].β*vr[lines[j].tail] )
isqold[j]=iold[j,1]^2 + iold[j,2]^2
#con=@constraint(m,isqnew[j]==2*(iold[j,1])*irij[j]+2*(iold[j,2])*ijij[j]-isqold[j])
if lines[j].Imax > 0
@constraint(m,-lines[j].Imax <= (irij[j]) <= lines[j].Imax)
@constraint(m,-lines[j].Imax <= (ijij[j]) <= lines[j].Imax)
end
#Infeasibility Handling for current flow from i to j
if iold[j,1]^2 + iold[j,2]^2 > lines[j].Imax^2
μi[j]=iold[j,2]/iold[j,1]
iold[j,1]=sign(iold[j,1])*(lines[j].Imax^2/(1+μi[j]^2))^0.5
iold[j,2]=sign(iold[j,2])*abs(μi[j] * iold[j,1])
@constraint(m,iold[j,1]* irij[j] + iold[j,2] * ijij[j] <= lines[j].Imax^2)
elseif iold[j,1]^2 + iold[j,2]^2 <= lines[j].Imax^2
@constraint(m,2*(iold[j,1])*irij[j]+2*(iold[j,2])*ijij[j]-isqold[j] <= lines[j].Imax^2)
end
#Constraint Reduction
#if (irijopt[h-1,j]^2+ijijopt[h-1,j]^2 >= r^2 * lines[j].Imax^2) | (irjiopt[h-1,j]^2+ijjiopt[h-1,j]^2 >= r^2 * lines[j].Imax^2)
end
for i = 1:numbuses
vsqold[i]=(vold[i,1])^2+(vold[i,2])^2
#@constraint(m,vsqnew[i]==2*(vold[i,1])*vr[i]+2*(vold[i,2])*vj[i]-vsqold[i])
#current injected at bus
@constraint(m,ir[i]== sum{irij[k], k in buses[i].outlist}-sum{irij[k], k in buses[i].inlist}+vr[i]*buses[i].Gs-vr[i]*buses[i].Bs)
@constraint(m, ij[i]==sum{ijij[k], k in buses[i].outlist}-sum{ijij[k], k in buses[i].inlist}+vj[i]*buses[i].Gs-vj[i]*buses[i].Bs)
@constraint(m,gp_tot[i]==sum{gp[k], k in buses[i].genids})
@constraint(m,gq_tot[i]==sum{gq[k], k in buses[i].genids})
#@constraint(m, sum{gp[k], k in buses[i].genids}-vold[i,1]*ir[i]-vold[i,2]*ij[i]-vr[i]*ibus_opt[i,1]-vj[i]*ibus_opt[i,2]+vold[i,1]*ibus_opt[i,1]+vold[i,2]*ibus_opt[i,2]==buses[i].Pd)
#@constraint(m, sum{gq[k], k in buses[i].genids}-vold[i,2]*ir[i]+vold[i,1]*ij[i]-vj[i]*ibus_opt[i,1]+vr[i]*ibus_opt[i,2]+vold[i,2]*ibus_opt[i,1]-vold[i,1]*ibus_opt[i,2]==buses[i].Qd)
@constraint(m, -buses[i].Vmax <= vr[i] <= buses[i].Vmax)
@constraint(m, -buses[i].Vmax <= vj[i] <= buses[i].Vmax)
#Infeasibility Handling
if vold[i,1]^2 + vold[i,2]^2 > buses[i].Vmax^2
μv[i]=vold[i,2]/vold[i,1]
vold[i,1]=sign(vold[i,1]*(buses[i].Vmax^2/(1+μv[i]^2))^0.5)
vold[i,2]=sign(vold[i,2])*abs(μv[i] * vold[i,1])
@constraint(m,vold[i,1]* vr[i] + vold[i,2] * vj[i] <= buses[i].Vmax^2)
elseif vold[i,1]^2 + vold[i,2]^2 <= buses[i].Vmax^2
@constraint(m,buses[i].Vmin^2 <= 2*(vold[i,1])*vr[i]+2*(vold[i,2])*vj[i]-vsqold[i] <= buses[i].Vmax^2)
end
P_bilinear[i] = vold[i,1]*iold[i,1]+vold[i,2]*iold[i,2]+buses[i].Pd
Q_bilinear[i] = vold[i,2]*iold[i,1]-vold[i,1]*iold[i,2]+buses[i].Qd
φpn[i] = max(0.001, abs(P_tay[i]-P_bilinear[i])/min(P_tay[i],P_bilinear[i]))
φqn[i] = max(0.001, abs(Q_tay[i]-Q_bilinear[i])/min(Q_tay[i],Q_bilinear[i]))
βr[i] = -0.15 * log10(φqn[i]) + 1
βj[i] = -0.15 * log10(φpn[i]) + 1
αr[i] = 1/βr[i]
αj[i] = 1/βj[i]
stepvr[i] = αr[i] * abs(buses[i].Vmax)/h^βr[i]
stepvj[i] = αj[i] * abs(buses[i].Vmax)/h^βj[i]
δp[i] = abs(P_tay[i]-P_bilinear[i])/min(P_tay[i],P_bilinear[i])
δq[i] = abs(Q_tay[i]-Q_bilinear[i])/min(Q_tay[i],Q_bilinear[i])
δp_tot=sum{δp[i],i=1:numbuses}
δq_tot=sum{δq[i],i=1:numbuses}
@constraint(m,-stepvr[i] <= vr[i]-vold[i,1] <= stepvr[i])
@constraint(m,-stepvj[i] <= vj[i]-vold[i,2] <= stepvj[i])
end
for i =1:numgenerators
@constraint(m,0<=gp[i]<=generators[i].Pgmax)
@constraint(m,-generators[i].Qgmax<=gq[i]<=generators[i].Qgmax)
end
solve(m)
return getvalue(gp),getvalue(gq),getobjectivevalue(m),getvalue(vr),getvalue(vj),getvalue(ir),getvalue(ij),getvalue(irij),getvalue(ijij),getvalue(gp_tot),getvalue(gq_tot)
end
(gp_opt,gq_opt,obj,vr_opt,vj_opt,ir_opt,ij_opt,irij_opt,ijij_opt,gp_tot_opt,gq_tot_opt)=solve_taylor(thermalLimitscale, generators, buses, lines);
println("**************************************************************************************")
println("The optimal cost is ", obj,"The iteration is", h,"vr is",vr_opt,"vj is",vj_opt)
println("**************************************************************************************")
#=vropt[h,:]= vr_opt
vjopt[h,:]= vj_opt
iropt[h,:]= ir_opt
ijopt[h,:]= ij_opt
irijopt[h,:]= irij_opt
ijijopt[h,:]= ijij_opt
=#
vold=hcat(vr_opt,vj_opt)
iold=hcat(irij_opt,ijij_opt)
ibus_opt=hcat(ir_opt,ij_opt)
P_tay=gp_tot_opt
Q_tay=gq_tot_opt
println(h)
h=h+1
end