lixiao00093
commented
1 month ago Is this a bug? I am quite sure this can work in ModelingToolkit v9.38. in the following version.
using ModelingToolkit, DifferentialEquations, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D_t
# --------system models --------
@connector AcBus begin
@variables begin
vm(t) = 1.0
va(t) = 0.0
pin(t) = 0.0, [connect = Flow]
qin(t) = 0.0, [connect = Flow]
end
end
@mtkmodel Zipload_component begin
@parameters begin
u = 1,[description = "Enable Signal"]
g
b
end
@components begin
bus = AcBus()
end
@equations begin
bus.pin ~ u * ( g * bus.vm ^2)
bus.qin ~ u * (- b * bus.vm ^2)
end
end
@mtkmodel Piline_component begin
@parameters begin
u = 1 ,[description = "Enable Signal"]
rft
xft
g0 = 0.0
b0
tap_ratio = 1.0
gft = rft/(rft^2 + xft^2)
bft = -xft/(rft^2 + xft^2)
m_real = tap_ratio
m2 = tap_ratio^2
end
@components begin
bus_fr = AcBus()
bus_to = AcBus()
end
@equations begin
bus_fr.pin*m2 ~ u*( g0/2*bus_fr.vm^2 + gft*bus_fr.vm^2
- gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
- bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
-bus_fr.qin*m2 ~ u*(b0/2*bus_fr.vm^2 + bft*bus_fr.vm^2
- bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
+ gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
bus_to.pin*m2 ~ u*( g0/2*m2*bus_to.vm^2 + gft*m2*bus_to.vm^2
- gft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
+ bft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
-bus_to.qin*m2 ~ u*(b0/2*m2*bus_to.vm^2 + bft*m2*bus_to.vm^2
- bft*m_real*bus_fr.vm*bus_to.vm*cos(bus_fr.va - bus_to.va)
- gft*m_real*bus_fr.vm*bus_to.vm*sin(bus_fr.va - bus_to.va))
end
end
@mtkmodel Syn2_component begin
@components begin
bus = AcBus()
end
@variables begin
delta(t), [description = "rotor angle"]
omega(t) = 1.0, [description = "rotor speed"]
vf(t)
tm(t)
p(t)
q(t)
_Id(t) = 0.0
_Iq(t) = 0.0
end
@parameters begin
time_flag
u = 1, [description = "Enable Signal"]
fn
H
Damp
ra
xd1
Kw
Kp
iM = u / (2 * H)
c1 = ra / (ra^2 + xd1^2)
c2 = xd1 / (ra^2 + xd1^2)
c3 = c2
vm0
va0
end
@equations begin
_Id ~ u*(-(bus.vm * sin(delta - bus.va) * c1) - (bus.vm * cos(delta - bus.va) * c3) + c3 * vf)
_Iq ~ u*( (bus.vm * sin(delta - bus.va) * c2) - (bus.vm * cos(delta - bus.va) * c1) + c1 * vf)
p ~ bus.vm * sin(delta - bus.va) * _Id + bus.vm * cos(delta - bus.va) * _Iq
q ~ bus.vm * cos(delta - bus.va) * _Id - bus.vm * sin(delta - bus.va) * _Iq
bus.pin ~ -p # for AcBus
bus.qin ~ -q # for AcBus
D_t(delta) ~ u * 2 * Ļ * fn * (omega - 1)
D_t(omega) ~ iM * (tm - p - ra * (_Id^2 + _Iq^2) - Damp * (omega - 1))
D_t(tm) + (time_flag<0)* (bus.va-va0) ~ 0
D_t(vf) + (time_flag<0)* (bus.vm-vm0) ~ 0
end
end
@mtkmodel Fault_component begin
@parameters begin
u = 0,[description="fault status"]
rf
xf
yg = rf/(rf^2 + xf^2)
yb = -xf/(rf^2 + xf^2)
end
@components begin
bus = AcBus()
end
@equations begin
bus.pin ~ u * yg * bus.vm^2
bus.qin ~ -u * yb * bus.vm^2
end
end
# -------colection of instances of models--------
Bus_Dict = Dict{Symbol,ODESystem}(Symbol(1)=>AcBus(;name = Symbol("Bus 1")),Symbol(2)=>AcBus(;name = Symbol("Bus 2")))
device_list::Dict{Symbol,ODESystem} = Dict(
:line => Piline_component(;name = :line, rft = 0.01,xft = 0.01,g0 = 0,b0 = 0.01,tap_ratio = 1.0),
:syn => Syn2_component(;name = :syn1, time_flag = -1, fn = 50.0, H = 1.0, Damp = 0.8, ra = 0.01, xd1 = 0.1, Kw=1.0, Kp=1.0, va0 = 0.0 , vm0=1.0),
:zipload => Zipload_component(;name = :zipload1, g = 1.0, b = 0.0),
:fault => Fault_component(;name = :fault, rf = 0.0, xf = 0.1)
)
# -------connection of models--------
eqs = [ connect(Bus_Dict[Symbol(1)], device_list[:syn].bus),
connect(Bus_Dict[Symbol(1)], device_list[:line].bus_fr),
connect(Bus_Dict[Symbol(1)], device_list[:fault].bus),
connect(Bus_Dict[Symbol(2)], device_list[:zipload].bus),
connect(Bus_Dict[Symbol(2)], device_list[:line].bus_to)]
# -------Apply and clear fault in the models--------
u_temp_record = Dict{Int,Vector{Float64}}()
function affect1!(integ, u, p, ctx)
if integ.ps[p.u] != ctx
println("apply fault at t = ", integ.t)
integ.ps[p.u] = ctx
u_temp_record[1] = copy(integ.u[integ.differential_vars.==false])
end
end
function affect2!(integ, u, p, ctx)
if integ.ps[p.u] != ctx
println("clear fault at t = ", integ.t)
integ.ps[p.u] = ctx
#offer the initial value of algebraic variables for the reinitialization
integ.u[integ.differential_vars.==false] = u_temp_record[1]
end
end
fault_function = [(t==1.0) => (affect1!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 1.0),
(t==1.05) => (affect2!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 0.0)]
# -------set the initial guessvalues of models--------
function set_ini_syn2(; ra, xd1)
pg0 = 1.0
qg0 = 0.0
vm0 = 1.0
va0 = 0.0
V = vm0 .* exp.(va0 .* 1im)
S = pg0 - qg0 .* 1im
I = S ./ conj.(V)
E = V + (ra + xd1 .* 1im) .* I
delta = imag(log.(E ./ abs.(E) + 0im))
delta0 = delta
omega0 = 1.0
# d- and q-axis voltages and currents
vdq = V .* 1im .* exp.(-delta .* 1im)
idq = I .* 1im .* exp.(-delta .* 1im)
vd = real(vdq)
vq = imag(vdq)
Id = real(idq)
Iq = imag(idq)
# mechanical torques/powers
tm0 = (vq + ra .* Iq) .* Iq + (vd + ra .* Id) .* Id
vf0 = vq + ra .* Iq + xd1 .* Id
return delta0, omega0, tm0, vf0, pg0, qg0
end
u0_guess = Dict(zip([device_list[:syn].delta, device_list[:syn].omega, device_list[:syn].tm, device_list[:syn].vf, device_list[:syn].p, device_list[:syn].q], set_ini_syn2(ra = 0.01, xd1 = 0.1)))
# -------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = fault_function)
test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)
test_system_xy0 = SteadyStateProblem(test_system_init_simplify, u0_guess; jac = true, sparse = true)
sol = DifferentialEquations.solve(test_system_xy0)
xy0_sol_dict = Dict{Symbolics.Num, Any}()
for var in unknowns(test_system_xy0.f.sys)
xy0_sol_dict[var] = sol[var]
end
for var in observed(test_system_xy0.f.sys)
xy0_sol_dict[var.lhs] = sol[var.lhs]
end
defaults(test_system_init_simplify)[device_list[:syn].time_flag] = 1
# -------format the system to solve--------
prob = ODEProblem(test_system_init_simplify, xy0_sol_dict, (0,20); fully_determined = true, jac=true,sparse=true)
# sol = solve(prob, Rodas4(), tstops = [1.0,1.05])
sol = solve(prob, dtmax = 0.01, tstops = [1.0,1.05])
plot(sol)
plot(sol.t,sol[device_list[:syn].delta])
plot(sol.t,sol[device_list[:syn].omega])
plot(sol.t,sol[device_list[:zipload].bus.vm])But this way cannot work in the later version, because the way of initialisation has changed. But I do not know why I can not run my code after I modify that like the top side of this issue.
When I run the above codes, the output of Pkg.status()
[6e4b80f9] BenchmarkTools v1.5.0
[336ed68f] CSV v0.10.14
ā [a93c6f00] DataFrames v1.6.1
[0c46a032] DifferentialEquations v7.14.0
[5789e2e9] FileIO v1.16.3
ā [961ee093] ModelingToolkit v9.38.0
[91a5bcdd] Plots v1.40.8
ā [0bca4576] SciMLBase v2.53.1
[37e2e46d] LinearAlgebra
[56ddb016] Logging
ChrisRackauckas
I think this relates to my previous issues: https://github.com/SciML/ModelingToolkit.jl/issues/3043, and has some relationship to https://github.com/SciML/OrdinaryDiffEq.jl/issues/2441.
I am stuck with this problem for a long time, no matter I use the latest version 9.41.0 or the older version of the ModelingToolkit.
What should I do? Should I modify my codes? Or this is a bug?
--------system models --------
@connector AcBus begin @variables begin vm(t), [guess = 1.0] va(t), [guess = 0.0] pin(t), [connect = Flow, guess = 0.0] qin(t), [connect = Flow, guess = 0.0] end end
@mtkmodel Zipload_component begin @parameters begin u = 1,[description = "Enable Signal"] g b end @components begin bus = AcBus() end @equations begin bus.pin ~ u ( g bus.vm ^2) bus.qin ~ u (- b bus.vm ^2) end end
@mtkmodel Piline_component begin @parameters begin u = 1 ,[description = "Enable Signal"] rft xft g0 = 0.0 b0 tap_ratio = 1.0 gft = rft/(rft^2 + xft^2) bft = -xft/(rft^2 + xft^2) m_real = tap_ratio m2 = tap_ratio^2 end
end
@mtkmodel Syn2_component begin @components begin bus = AcBus() end
end
@mtkmodel Fault_component begin @parameters begin u = 0,[description="fault status"] rf xf yg = rf/(rf^2 + xf^2) yb = -xf/(rf^2 + xf^2) end
end
-------colection of instances of models--------
Bus_Dict = Dict{Symbol,ODESystem}(Symbol(1)=>AcBus(;name = Symbol("Bus 1")),Symbol(2)=>AcBus(;name = Symbol("Bus 2"))) device_list::Dict{Symbol,ODESystem} = Dict( :line => Piline_component(;name = :line, rft = 0.01,xft = 0.01,g0 = 0,b0 = 0.01,tap_ratio = 1.0), :syn => Syn2_component(;name = :syn1, fn = 50.0, H = 1.0, Damp = 0.8, ra = 0.01, xd1 = 0.1, Kw=1.0, Kp=1.0, va0 = 0.0 , vm0=1.0), :zipload => Zipload_component(;name = :zipload1, g = 1.0, b = 0.0), :fault => Fault_component(;name = :fault, rf = 0.0, xf = 0.1) )
-------connection of models--------
eqs = [ connect(Bus_Dict[Symbol(1)], device_list[:syn].bus), connect(Bus_Dict[Symbol(1)], device_list[:line].bus_fr), connect(Bus_Dict[Symbol(1)], device_list[:fault].bus),
-------initialization of models--------
ini_eqs = [device_list[:syn].omega ~ 1, device_list[:syn].tm - device_list[:syn].p - device_list[:syn].ra * (device_list[:syn]._Id^2 + device_list[:syn]._Iq^2) ~ 0, device_list[:syn].bus.va ~ device_list[:syn].va0, device_list[:syn].bus.vm ~ device_list[:syn].vm0]
-------change the parameters of models--------
-------Apply and clear fault in the models--------
u_temp_record = Dict{Int,Vector{Float64}}() function affect1!(integ, u, p, ctx) println("apply fault at t = ", integ.t) if integ.ps[p.u] != ctx integ.ps[p.u] = ctx u_temp_record[1] = copy(integ.u[integ.differential_vars.==false]) end end function affect2!(integ, u, p, ctx) println("clear fault at t = ", integ.t) if integ.ps[p.u] != ctx integ.ps[p.u] = ctx
offer the initial value of algebraic variables for the reinitialization
end
fault_function = [(t==1.0) => (affect1!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 1.0), (t==1.05) => (affect2!, [], [device_list[:fault].u=>:u], [device_list[:fault].u], 0.0)]
-------set the initial guessvalues of models--------
function set_ini_syn2(; ra, xd1) pg0 = 1.0 qg0 = 0.0 vm0 = 1.0 va0 = 0.0
end u0_guess = Dict(zip([device_list[:syn].delta, device_list[:syn].omega, device_list[:syn].tm, device_list[:syn].vf, device_list[:syn].p, device_list[:syn].q], set_ini_syn2(ra = 0.01, xd1 = 0.1)))
-------simplify the system--------
test_system_init = ODESystem(eqs, t,[],[]; systems = collect(values(merge(device_list, Bus_Dict))), name = :test_system_init, discrete_events = fault_function) test_system_init_simplify = structural_simplify(test_system_init, fully_determined = true)
-------format the system to solve--------
prob = ODEProblem(test_system_init_simplify,[], (0,20),[]; guesses = u0_guess, initialization_eqs = ini_eqs, fully_determined = true, jac=true,sparse=true)
sol = solve(prob, Rodas4(), tstops = [1.0,1.05])
sol = solve(prob, tstops = [1.0,1.05]) plot(sol)
[6e4b80f9] BenchmarkTools v1.5.0 [336ed68f] CSV v0.10.14 [a93c6f00] DataFrames v1.7.0 [0c46a032] DifferentialEquations v7.14.0 [5789e2e9] FileIO v1.16.3 [961ee093] ModelingToolkit v9.41.0 [8913a72c] NonlinearSolve v3.14.0 [1dea7af3] OrdinaryDiffEq v6.89.0 [91a5bcdd] Plots v1.40.8 [0bca4576] SciMLBase v2.54.0 [37e2e46d] LinearAlgebra [56ddb016] Logging [2f01184e] SparseArrays v1.10.0
The output:
if solving the problem with discrete event by Rodas4() solver:
The output: