Dane Sabo
ed59b30dd0
D Presentations/ERLM/bouncing_ball_hybrid.png A Presentations/ERLM/images/4_research_approach/phase_portrait_sg1.png A Presentations/ERLM/images/4_research_approach/two_loop.png A Presentations/ERLM/images/7_broader_impacts/billion.jpg M Presentations/ERLM/main.aux M Presentations/ERLM/main.fdb_latexmk M Presentations/ERLM/main.fls M Presentations/ERLM/main.log
216 lines
6.9 KiB
Julia
216 lines
6.9 KiB
Julia
using ReachabilityAnalysis
|
|
using LinearAlgebra
|
|
using Plots
|
|
|
|
"""
|
|
Two-Loop Reactor: Linearized Reachability Analysis
|
|
|
|
Since the nonlinear system is causing issues with all algorithms,
|
|
let's linearize around an equilibrium point and use LGG09.
|
|
|
|
This is actual formal reachability analysis - the linearization is
|
|
a standard approximation technique used in control theory.
|
|
"""
|
|
|
|
println("=== Two-Loop Reactor: Linearized Reachability ===\n")
|
|
|
|
# Constants
|
|
const C_0 = 33.33
|
|
const τ_0 = 0.75 * C_0
|
|
const μ = 0.6
|
|
const P_r = 100.0
|
|
const Q2 = 50.0
|
|
|
|
C_r = μ * C_0
|
|
C_sg = (1 - μ) * C_0 / 2
|
|
W = C_0 / (2 * τ_0)
|
|
|
|
println("=== System Parameters ===")
|
|
println("C_r = $C_r, C_sg = $C_sg, W = $W")
|
|
|
|
# Find equilibrium with Q1 = 50
|
|
# At equilibrium: dT/dt = 0
|
|
# From equations:
|
|
# P_r - W*(T_hot - T_cold1) - W*(T_hot - T_cold2) = 0
|
|
# W*(T_hot - T_cold1) - Q1 = 0 => T_hot - T_cold1 = Q1/W
|
|
# W*(T_hot - T_cold2) - Q2 = 0 => T_hot - T_cold2 = Q2/W
|
|
|
|
Q1_nominal = 50.0
|
|
T_hot_eq = 453.0 # Approximate equilibrium
|
|
T_cold1_eq = T_hot_eq - Q1_nominal/W
|
|
T_cold2_eq = T_hot_eq - Q2/W
|
|
|
|
println("\n=== Equilibrium Point (Q1 = $Q1_nominal) ===")
|
|
println("T_hot_eq = $T_hot_eq °F")
|
|
println("T_cold1_eq = $T_cold1_eq °F")
|
|
println("T_cold2_eq = $T_cold2_eq °F")
|
|
|
|
# Linearize around equilibrium
|
|
# State: x = [T_hot - T_hot_eq, T_cold1 - T_cold1_eq, T_cold2 - T_cold2_eq, Q1 - Q1_nominal]
|
|
#
|
|
# Original equations:
|
|
# dT_hot/dt = (P_r - W*(T_hot - T_cold1) - W*(T_hot - T_cold2)) / C_r
|
|
# dT_cold1/dt = (W*(T_hot - T_cold1) - Q1) / C_sg
|
|
# dT_cold2/dt = (W*(T_hot - T_cold2) - Q2) / C_sg
|
|
# dQ1/dt = 0
|
|
#
|
|
# Linearized (using deviations from equilibrium):
|
|
# dx1/dt = (-W*x1 + W*x2 + W*x3) / C_r
|
|
# dx2/dt = (W*x1 - W*x2 - x4) / C_sg
|
|
# dx3/dt = (W*x1 - W*x3) / C_sg
|
|
# dx4/dt = 0
|
|
|
|
A = [
|
|
-W/C_r W/C_r W/C_r 0.0;
|
|
W/C_sg -W/C_sg 0.0 -1.0/C_sg;
|
|
W/C_sg 0.0 -W/C_sg 0.0;
|
|
0.0 0.0 0.0 0.0
|
|
]
|
|
|
|
println("\n=== Linearized System Matrix A ===")
|
|
display(A)
|
|
println()
|
|
|
|
# Initial set: Small deviations from equilibrium
|
|
# Starting close to equilibrium but with full Q1 uncertainty
|
|
δT_init = 2.0 # ±2°F initial temperature deviation
|
|
Q1_min = 45.0
|
|
Q1_max = 55.0
|
|
|
|
X0 = Hyperrectangle(
|
|
low=[-δT_init, -δT_init, -δT_init, Q1_min - Q1_nominal],
|
|
high=[δT_init, δT_init, δT_init, Q1_max - Q1_nominal]
|
|
)
|
|
|
|
println("\n=== Initial Set (deviations from equilibrium) ===")
|
|
println("δT_hot ∈ [-$δT_init, $δT_init]")
|
|
println("δT_cold1 ∈ [-$δT_init, $δT_init]")
|
|
println("δT_cold2 ∈ [-$δT_init, $δT_init]")
|
|
println("δQ1 ∈ [$(Q1_min - Q1_nominal), $(Q1_max - Q1_nominal)]")
|
|
|
|
# Create linear system
|
|
sys = @system(x' = A*x, x ∈ Universe(4))
|
|
|
|
# Create IVP
|
|
prob = InitialValueProblem(sys, X0)
|
|
|
|
println("\n=== Solving with LGG09 (Exact for Linear Systems!) ===")
|
|
|
|
# Solve - LGG09 is EXACT for linear systems!
|
|
# Must specify vars and dimension
|
|
sol = solve(prob, T=30.0, alg=LGG09(δ=0.5, vars=[1, 2, 3, 4], dim=4))
|
|
|
|
println("✓ Success! Computed $(length(sol)) reach sets")
|
|
println("\nNote: LGG09 provides EXACT reachability for linear systems!")
|
|
|
|
# Convert back to absolute temperatures for plotting
|
|
println("\n=== Creating Plots ===")
|
|
|
|
# Extract and shift back to absolute coordinates
|
|
function plot_absolute(sol, var_idx, eq_val, ylabel_text, title_text, color_choice)
|
|
# Get the reach sets
|
|
times = Float64[]
|
|
lower = Float64[]
|
|
upper = Float64[]
|
|
|
|
for reach_set in sol
|
|
t = inf(reach_set.t_start)
|
|
set = reach_set.X
|
|
|
|
# Project to variable
|
|
proj = project(set, var_idx)
|
|
|
|
# Get bounds
|
|
low_val = low(proj)[1] + eq_val
|
|
high_val = high(proj)[1] + eq_val
|
|
|
|
push!(times, t)
|
|
push!(lower, low_val)
|
|
push!(upper, high_val)
|
|
end
|
|
|
|
p = plot(times, lower, fillrange=upper,
|
|
xlabel="Time (s)", ylabel=ylabel_text,
|
|
title=title_text,
|
|
fillalpha=0.5, color=color_choice, lw=2,
|
|
label="Reach tube", legend=:topright)
|
|
return p
|
|
end
|
|
|
|
p1 = plot_absolute(sol, 1, T_hot_eq, "T_hot (°F)", "Hot Leg Temperature", :red)
|
|
p2 = plot_absolute(sol, 2, T_cold1_eq, "T_cold1 (°F)", "Cold Leg 1 (SG1 - Uncertain)", :blue)
|
|
p3 = plot_absolute(sol, 3, T_cold2_eq, "T_cold2 (°F)", "Cold Leg 2 (SG2 - Known)", :green)
|
|
|
|
# Phase portrait: T_hot vs T_cold1
|
|
p4 = plot(xlabel="T_cold1 (°F)", ylabel="T_hot (°F)",
|
|
title="Phase Portrait: T_hot vs T_cold1",
|
|
legend=:topright)
|
|
|
|
for reach_set in sol
|
|
set = reach_set.X
|
|
|
|
# Project to (x2, x1)
|
|
proj = project(set, [2, 1])
|
|
|
|
# Get vertices and shift
|
|
verts = vertices_list(proj)
|
|
verts_abs = [[v[1] + T_cold1_eq, v[2] + T_hot_eq] for v in verts]
|
|
|
|
# Plot polygon
|
|
xs = [v[1] for v in verts_abs]
|
|
ys = [v[2] for v in verts_abs]
|
|
push!(xs, xs[1]) # Close the polygon
|
|
push!(ys, ys[1])
|
|
|
|
plot!(p4, xs, ys, fillalpha=0.3, color=:purple, lw=0, label="")
|
|
end
|
|
plot!(p4, [], [], fillalpha=0.3, color=:purple, label="Reachable sets")
|
|
|
|
plot_combined = plot(p1, p2, p3, p4, layout=(2, 2), size=(1400, 1000),
|
|
plot_title="Linearized Reachability Analysis: Q1 ∈ [$Q1_min, $Q1_max]")
|
|
savefig(plot_combined, "two_loop_linearized_reachability.png")
|
|
println("Saved: two_loop_linearized_reachability.png")
|
|
|
|
# Detailed phase portrait
|
|
p_phase = plot(xlabel="T_cold1 (°F)", ylabel="T_hot (°F)",
|
|
title="Phase Portrait: T_hot vs T_cold1\n(Linearized Reachability: Q1 ∈ [$Q1_min, $Q1_max])",
|
|
size=(1000, 900), legend=:topright)
|
|
|
|
for (i, reach_set) in enumerate(sol)
|
|
set = reach_set.X
|
|
proj = project(set, [2, 1])
|
|
verts = vertices_list(proj)
|
|
verts_abs = [[v[1] + T_cold1_eq, v[2] + T_hot_eq] for v in verts]
|
|
|
|
xs = [v[1] for v in verts_abs]
|
|
ys = [v[2] for v in verts_abs]
|
|
push!(xs, xs[1])
|
|
push!(ys, ys[1])
|
|
|
|
alpha_val = 0.2 + 0.3 * (i / length(sol)) # Fade over time
|
|
plot!(p_phase, xs, ys, fillalpha=alpha_val, color=:purple, lw=0.5, label="")
|
|
end
|
|
|
|
plot!(p_phase, [], [], fillalpha=0.5, color=:purple, label="Reachable region")
|
|
|
|
savefig(p_phase, "phase_portrait_linearized.png")
|
|
println("Saved: phase_portrait_linearized.png")
|
|
|
|
println("\n✓ Complete!")
|
|
println("\n=== Results ===")
|
|
println("This is FORMAL reachability analysis using linear approximation.")
|
|
println("LGG09 provides EXACT reachable sets for linear systems.")
|
|
println("\nThe linearization is valid near the equilibrium point.")
|
|
println("For large deviations, use parameter sweeps (our previous approach).")
|
|
|
|
# Print final bounds
|
|
final_set = sol[end].X
|
|
println("\n=== Final Reachable Set (t=30s) ===")
|
|
proj1 = project(final_set, 1)
|
|
proj2 = project(final_set, 2)
|
|
proj3 = project(final_set, 3)
|
|
|
|
println("T_hot ∈ [$(round(low(proj1)[1] + T_hot_eq, digits=1)), $(round(high(proj1)[1] + T_hot_eq, digits=1))] °F")
|
|
println("T_cold1 ∈ [$(round(low(proj2)[1] + T_cold1_eq, digits=1)), $(round(high(proj2)[1] + T_cold1_eq, digits=1))] °F")
|
|
println("T_cold2 ∈ [$(round(low(proj3)[1] + T_cold2_eq, digits=1)), $(round(high(proj3)[1] + T_cold2_eq, digits=1))] °F")
|