Open
Description
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)?