ModelingToolkit.jl
ModelingToolkit.jl copied to clipboard
Transform NonlinearSystem to OptimizationSystem?
I have a quation about transform NonlinearSystem and OptimizationSystem. I construct a system based on NonlinearSystem,a source and two turbine and a sink, change source.p and sink.p,we can get the flow through the system. if I want fix the flow to some const value, modify the source.p or sink.p, then need construct the optimization problem. so I want to know can I transform the nonlinearsystem to optimizationsystem, the I provide the loss function like loss = abs(turbine1.a.m - 700.)
, and solve the optimization problem in ModelingToolkit.
using XSteam
using ModelingToolkit
using NonlinearSolve
using GalacticOptim
@register XSteam.phs(p, h)
@register XSteam.psh(p, s)
@register XSteam.pht(p, h)
@register XSteam.phs(p, h)
@register XSteam.phv(p, h)
function Port(;name)
@variables p m h
NonlinearSystem(Equation[], [p, m, h], [], name=name, defaults=Dict(p=>1.0, m=>0.0, h=>1000))
end
function Source(;name, p, h)
val_p = p
val_h = h
@named a = Port()
@parameters p h
eqs = [
a.p ~ p
a.h ~ h
]
NonlinearSystem(eqs, [], [p, h], systems=[a], defaults=Dict(p => val_p, h => val_h), name=name)
end
function Sink(;name, p)
val_p = p
@named a = Port()
@parameters p
eqs = [
a.p ~ p
]
NonlinearSystem(eqs, [], [p], systems=[a], defaults=Dict(p => val_p), name=name)
end
function connect_ports(ps...)
eqs = [
0 ~ sum(p->p.m, ps) # sum(mdot) = 0
]
for i in 2:length(ps)
push!(eqs, ps[1].p ~ ps[i].p)
push!(eqs, ps[1].h ~ ps[i].h)
end
return eqs
end
function Turbine(;name, p_in_ref, p_out_ref, h_in_ref, η_ref, m_ref)
val_p₀ = p_in_ref
val_p₁ = p_out_ref
val_h₀ = h_in_ref
val_η = η_ref
val_m₀ = m_ref
@named a = Port()
@named b = Port()
@variables power
@parameters p_in_ref,p_out_ref,h_in_ref,η_ref,m_ref
eqs = [
0 ~ a.m - m_ref*sqrt((a.p^2 - b.p^2)/(p_in_ref^2 - p_out_ref^2)) * sqrt(p_in_ref * XSteam.phv(a.p, a.h) / a.p / XSteam.phv(p_in_ref, h_in_ref))
η ~ (a.h - b.h) / (a.h - XSteam.psh(b.p, XSteam.phs(a.p, a.h)))
0 ~ a.m + b.m
0 ~ power - a.m * (a.h - b.h)
]
NonlinearSystem(eqs,[power], [p_in_ref, p_out_ref, h_in_ref, η_ref, m_ref], systems=[a, b]
, defaults=Dict(power => 0.0, p_in_ref => val_p₀, p_out_ref => val_p₁, h_in_ref => val_h₀, η_ref => val_η, m_ref => val_m₀), name=name)
end
const p₀ = 270
const h₀ = 3475.1
const p₁ = 79.09
const m₀ = 800.253
const η = 0.8831
@named source = Source(p = 200.0, h = h₀)
@named turbine1 = Turbine(p_in_ref = p₀, p_out_ref = p₁, h_in_ref = h₀, η_ref = η, m_ref = m₀)
@named turbine2 = Turbine(p_in_ref = 79.09, p_out_ref = 60.0, h_in_ref = 3069.3, η_ref = 0.918, m_ref = m₀ -41.0)
@named sink = Sink(p = 60.0)
cnnt_eqs = [
connect_ports(source.a, turbine1.a)
connect_ports(turbine1.b, turbine2.a)
connect_ports(turbine2.b, sink.a)
]
@named tes_model = NonlinearSystem(cnnt_eqs, [], [], systems=[source, turbine1, turbine2, sink])
sys = structural_simplify(tes_model)
u0 = [
turbine1.b.p => 65.0
turbine1.b.h => 3106.0
turbine1.b.h => 3039.0
]
# these some way to transfrom the sys to a OptimizationSystem?
prob = NonlinearProblem(sys, u0)
sol = solve(prob, NewtonRaphson())
Did OptimizationSystem have signature like OptimizatonSystem(sys,loss,u,p)?
What would this system do? Have the nonlinear systems be the constraints?
the system calculate the flow throught the turbine according the input and output pressure. yes, the nonlinear system be the optimization system constraint, can implement in ModelingToolkit?
Not right now.