From 018f7804e1bc5eac7a2d0479bbf7d3e15c42341f Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Wed, 4 Sep 2024 18:15:55 -0700
Subject: [PATCH 1/6] Risk expressions + objectives

---
 src/constants.jl           |  1 +
 src/expressions/dynamic.jl | 12 ++++++++++++
 src/extract.jl             | 10 +++++-----
 src/objectives.jl          |  2 ++
 4 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/constants.jl b/src/constants.jl
index 8407c6f4..be328604 100644
--- a/src/constants.jl
+++ b/src/constants.jl
@@ -134,6 +134,7 @@ const index_2_name__stability__model = Dict(
     202 => :q08_gt_2, # q(@rho=0.8) > 2.
     # 300s: Density Limit Models
     301 => :gw_density, # Density limit defined by Greenwald fraction
+    302 => :edge_collisionality,
     # 400s: Shaping Limit Models
     401 => :κ_controllability,
     # 900s: Stability Codes
diff --git a/src/expressions/dynamic.jl b/src/expressions/dynamic.jl
index d6c4ff33..5052b9e3 100644
--- a/src/expressions/dynamic.jl
+++ b/src/expressions/dynamic.jl
@@ -590,6 +590,18 @@ dyexp["pulse_schedule.tf.b_field_tor_vacuum.reference"] =
 dyexp["pulse_schedule.tf.time"] =
     (time; dd, _...) ->  dd.equilibrium.time
 
+#= ==== =#
+#  risk  #
+#= ==== =#
+dyexp["risk.engineering.loss[:].risk"] = 
+    (; loss, _...) -> loss.probability * loss.severity 
+
+dyexp["risk.plasma.risk"] = 
+    (; dd, _...) -> isempty(dd.risk.plasma) ? 0.0 : sum(loss.risk for loss in dd.risk.plasma.loss)    
+
+dyexp["risk.engineering.risk"] = 
+    (; dd, _...) -> isempty(dd.risk.engineering) ? 0.0 : sum(loss.risk for loss in dd.risk.engineering.loss)
+
 #= ========= =#
 #  stability  #
 #= ========= =#
diff --git a/src/extract.jl b/src/extract.jl
index ca92126d..0300b36d 100644
--- a/src/extract.jl
+++ b/src/extract.jl
@@ -127,10 +127,10 @@ function update_ExtractFunctionsLibrary!()
     ExtractLibFunction(:build, :OH_stress_margin, "-", dd -> dd.solid_mechanics.center_stack.properties.yield_strength.oh/maximum(dd.solid_mechanics.center_stack.stress.vonmises.oh))
 
     ExtractLibFunction(:costing, :levelized_CoE, "\$/kWh", dd -> dd.costing.levelized_CoE)
-    ExtractLibFunction(:costing, :TF_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"TF") / dd.costing.cost_direct_capital.cost)
-    ExtractLibFunction(:costing, :BOP_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"balance of plant equipment") / dd.costing.cost_direct_capital.cost)
-    ExtractLibFunction(:costing, :blanket_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"blanket") / dd.costing.cost_direct_capital.cost)
-    ExtractLibFunction(:costing, :cryostat_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"cryostat") / dd.costing.cost_direct_capital.cost)
+    ExtractLibFunction(:costing, :TF_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"TF") / dd.costing.cost_direct_capital.cost)
+    ExtractLibFunction(:costing, :BOP_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"balance of plant equipment") / dd.costing.cost_direct_capital.cost)
+    ExtractLibFunction(:costing, :blanket_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"blanket") / dd.costing.cost_direct_capital.cost)
+    ExtractLibFunction(:costing, :cryostat_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"cryostat") / dd.costing.cost_direct_capital.cost)
 
     ExtractLibFunction(:costing, :capital_cost, "\$B", dd -> dd.costing.cost_direct_capital.cost / 1E3)
     ExtractLibFunction(:constraint, :min_required_power_electric_net, "-", dd -> CFL[:min_required_power_electric_net](dd))
@@ -308,7 +308,7 @@ function print_tiled(io::IO, xtract::AbstractDict{Symbol,ExtractFunction}; termi
     end
 end
 
-function select_direct_captial_cost(dd::IMAS.dd, what::String)
+function select_direct_capital_cost(dd::IMAS.dd, what::String)
     for sys in dd.costing.cost_direct_capital.system
         idx = findfirst(x-> x.name ==what, sys.subsystem)
         if !isnothing(idx)
diff --git a/src/objectives.jl b/src/objectives.jl
index c57687ce..126a03c9 100644
--- a/src/objectives.jl
+++ b/src/objectives.jl
@@ -27,6 +27,8 @@ function update_ObjectiveFunctionsLibrary!()
     ObjectiveFunction(:max_log10_flattop, "log₁₀(hours)", dd -> log10(dd.build.oh.flattop_duration / 3600.0), Inf)
     ObjectiveFunction(:min_βn, "", dd -> dd.equilibrium.time_slice[].global_quantities.beta_normal, -Inf)
     ObjectiveFunction(:min_R0, "m", dd -> dd.equilibrium.time_slice[].boundary.geometric_axis.r, -Inf)
+    ObjectiveFunction(:min_engineering_risk, "", dd -> dd.risk.engineering.risk, -Inf)
+    ObjectiveFunction(:min_plasma_risk, "", dd -> dd.risk.plasma.risk, -Inf)
     #! format: on
     return ObjectiveFunctionsLibrary
 end

From 57516379a588d3bc70253dd486ebce1511fc2eae Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Thu, 12 Sep 2024 10:48:03 -0700
Subject: [PATCH 2/6] Extract relevant risk results

---
 src/extract.jl    | 3 +++
 src/objectives.jl | 4 ++--
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/extract.jl b/src/extract.jl
index 6ad115cf..b46325eb 100644
--- a/src/extract.jl
+++ b/src/extract.jl
@@ -135,6 +135,9 @@ function update_ExtractFunctionsLibrary!()
     ExtractLibFunction(:costing, :blanket_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"blanket") / dd.costing.cost_direct_capital.cost)
     ExtractLibFunction(:costing, :cryostat_of_total, "%", dd -> 100 * select_direct_capital_cost(dd,"cryostat") / dd.costing.cost_direct_capital.cost)
 
+    ExtractLibFunction(:risk, :total_engineering_risk, "\$M", dd -> dd.risk.engineering.risk)
+    ExtractLibFunction(:risk, :total_plasma_risk, "\$/kWh", dd -> dd.risk.plasma.risk)
+
     ExtractLibFunction(:costing, :capital_cost, "\$B", dd -> dd.costing.cost_direct_capital.cost / 1E3)
     ExtractLibFunction(:constraint, :min_required_power_electric_net, "-", dd -> CFL[:min_required_power_electric_net](dd))
     ExtractLibFunction(:constraint, :required_power_electric_net, "-", dd -> CFL[:required_power_electric_net](dd))
diff --git a/src/objectives.jl b/src/objectives.jl
index 126a03c9..0ecb6999 100644
--- a/src/objectives.jl
+++ b/src/objectives.jl
@@ -27,8 +27,8 @@ function update_ObjectiveFunctionsLibrary!()
     ObjectiveFunction(:max_log10_flattop, "log₁₀(hours)", dd -> log10(dd.build.oh.flattop_duration / 3600.0), Inf)
     ObjectiveFunction(:min_βn, "", dd -> dd.equilibrium.time_slice[].global_quantities.beta_normal, -Inf)
     ObjectiveFunction(:min_R0, "m", dd -> dd.equilibrium.time_slice[].boundary.geometric_axis.r, -Inf)
-    ObjectiveFunction(:min_engineering_risk, "", dd -> dd.risk.engineering.risk, -Inf)
-    ObjectiveFunction(:min_plasma_risk, "", dd -> dd.risk.plasma.risk, -Inf)
+    ObjectiveFunction(:min_engineering_risk, "\$M", dd -> dd.risk.engineering.risk, -Inf)
+    ObjectiveFunction(:min_plasma_risk, "\$/kWh", dd -> dd.risk.plasma.risk, -Inf)
     #! format: on
     return ObjectiveFunctionsLibrary
 end

From 012d54d3f205e252741b52c5cbf1ceff19dea125 Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Thu, 3 Oct 2024 12:16:51 -0700
Subject: [PATCH 3/6] Basic risk plots

---
 src/plot.jl | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)

diff --git a/src/plot.jl b/src/plot.jl
index bf05f616..0f537e43 100644
--- a/src/plot.jl
+++ b/src/plot.jl
@@ -407,6 +407,67 @@ end
     end
 end
 
+# ==== #
+# risk #
+# ==== #
+
+@recipe function plot_risk(rsk::IMAS.risk)
+    eng_loss = rsk.engineering.loss
+    eng_names = ["$(sys_loss.description)" for sys_loss in reverse(eng_loss)]
+    eng_risks = [sys_loss.risk for sys_loss in reverse(eng_loss)]
+    eng_perc = ["$(round(sys_loss.risk/sum(eng_risks)*100))%" for sys_loss in reverse(eng_loss)]
+
+    plasma_loss = rsk.plasma.loss
+    plasma_names = ["$(sys_loss.description)" for sys_loss in reverse(plasma_loss)]
+    plasma_risks = [sys_loss.risk for sys_loss in reverse(plasma_loss)]
+    plasma_perc = ["$(round(sys_loss.risk/sum(plasma_risks)*100))%" for sys_loss in reverse(plasma_loss)]
+
+    size --> (800, 400)
+    layout := RecipesBase.@layout (1, 2)
+
+    @series begin 
+        subplot := 1
+        seriestype := :bar
+        orientation := :horizontal
+        title := "Engineering risk" * "   " * @sprintf("[Total = %.3g \$M]", sum(eng_risks))
+        titlefontsize := 10
+        ylim := (0, length(eng_risks))
+        label := ""
+        annotation := [(0.0, kk - 0.5, ("  $x  $(titlecase(n,strict=false))", :left, 8)) for (kk, (c, x, n)) in enumerate((collect(zip(eng_risks, eng_perc, eng_names))))]
+        annotationvalign := :center
+        label := ""
+        xticks := 0:round(maximum(eng_risks) / 4, digits = 1):round(maximum(eng_risks), digits = 1)
+        xlabel := "[\$M]"
+        showaxis := :x
+        yaxis := nothing
+        alpha := 0.5
+        linecolor := :match
+        color := PlotUtils.palette(:tab10)[1]
+        eng_names, eng_risks
+    end
+
+    @series begin 
+        subplot := 2
+        seriestype := :bar
+        orientation := :horizontal
+        title := "Plasma risk" * "    " * @sprintf("[Total = %.3g \$/kWh]", sum(plasma_risks))
+        titlefontsize := 10
+        ylim := (0, length(plasma_risks))
+        label := ""
+        annotation := [(0.0, kk - 0.5, ("  $x  $(titlecase(n,strict=false))", :left, 8)) for (kk, (c, x, n)) in enumerate((collect(zip(plasma_risks, plasma_perc, plasma_names))))]
+        annotationvalign := :center
+        label := ""
+        xticks := 0:round(maximum(plasma_risks) / 4, digits = 1):round(maximum(plasma_risks), digits = 1)
+        xlabel := "[\$/kWh]"
+        showaxis := :x
+        yaxis := nothing
+        alpha := 0.5
+        linecolor := :match
+        color := PlotUtils.palette(:tab10)[2]
+        plasma_names, plasma_risks
+    end
+end
+
 # =========== #
 # equilibrium #
 # =========== #

From 47cbee3d44770da87449e44e9156a9eb4d9ec049 Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Mon, 21 Oct 2024 13:56:47 -0700
Subject: [PATCH 4/6] Add plots with risks as error on cost

---
 src/plot.jl | 40 ++++++++++++++++++++++++++++++++++++++--
 1 file changed, 38 insertions(+), 2 deletions(-)

diff --git a/src/plot.jl b/src/plot.jl
index 0f537e43..77ba64ba 100644
--- a/src/plot.jl
+++ b/src/plot.jl
@@ -422,8 +422,12 @@ end
     plasma_risks = [sys_loss.risk for sys_loss in reverse(plasma_loss)]
     plasma_perc = ["$(round(sys_loss.risk/sum(plasma_risks)*100))%" for sys_loss in reverse(plasma_loss)]
 
-    size --> (800, 400)
-    layout := RecipesBase.@layout (1, 2)
+    dd = parent(rsk)
+    direct_capital_cost = dd.costing.cost_direct_capital.cost 
+    levelized_cost = dd.costing.levelized_CoE
+
+    size --> (800, 800)
+    layout := RecipesBase.@layout (2, 2)
 
     @series begin 
         subplot := 1
@@ -437,6 +441,7 @@ end
         annotationvalign := :center
         label := ""
         xticks := 0:round(maximum(eng_risks) / 4, digits = 1):round(maximum(eng_risks), digits = 1)
+        xlim := (0, maximum(eng_risks))
         xlabel := "[\$M]"
         showaxis := :x
         yaxis := nothing
@@ -458,6 +463,7 @@ end
         annotationvalign := :center
         label := ""
         xticks := 0:round(maximum(plasma_risks) / 4, digits = 1):round(maximum(plasma_risks), digits = 1)
+        xlim := (0, maximum(plasma_risks))
         xlabel := "[\$/kWh]"
         showaxis := :x
         yaxis := nothing
@@ -466,6 +472,36 @@ end
         color := PlotUtils.palette(:tab10)[2]
         plasma_names, plasma_risks
     end
+
+    @series begin 
+        subplot := 3
+        seriestype := :bar 
+        label := "Direct capital cost with error"
+        ylabel := "Direct capital cost (\$M)"
+        orientation := :vertical 
+        legend := :bottomleft
+        title := "Error on direct capital cost"
+        alpha := 0.5
+        yerror := [rsk.engineering.risk]
+        color := PlotUtils.palette(:tab10)[1]
+        [""], [direct_capital_cost]
+    end
+    
+    @series begin
+        subplot := 4
+        seriestype := :bar
+        label := "Levelized cost with error"
+        ylabel := "Levelized cost of electricity (\$/kWh)"
+        orientation := :vertical
+        legend := :bottomleft
+        title := "Error on levelized cost"
+        right_ylabel := "Levelized cost of electricity with error"
+        alpha := 0.5
+        yerror := [rsk.plasma.risk]
+        color := PlotUtils.palette(:tab10)[2]
+        [""], [levelized_cost]
+    end
+
 end
 
 # =========== #

From 910311c8032f365abb0114c44b7e3c5b86630988 Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Wed, 5 Mar 2025 14:13:57 -0700
Subject: [PATCH 5/6] Update expressions for failure modes in engineering risk

---
 src/expressions/dynamic.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/expressions/dynamic.jl b/src/expressions/dynamic.jl
index be9cd0b2..81610a98 100644
--- a/src/expressions/dynamic.jl
+++ b/src/expressions/dynamic.jl
@@ -615,7 +615,7 @@ dyexp["pulse_schedule.lh.power.reference"] =
 #  risk  #
 #= ==== =#
 dyexp["risk.engineering.loss[:].risk"] = 
-    (; loss, _...) -> loss.probability * loss.severity 
+    (; loss, _...) -> (sum((fm.probability * fm.weight) for fm in loss.failure_mode) * loss.severity)
 
 dyexp["risk.plasma.risk"] = 
     (; dd, _...) -> isempty(dd.risk.plasma) ? 0.0 : sum(loss.risk for loss in dd.risk.plasma.loss)    
@@ -623,7 +623,7 @@ dyexp["risk.plasma.risk"] =
 dyexp["risk.engineering.risk"] = 
     (; dd, _...) -> isempty(dd.risk.engineering) ? 0.0 : sum(loss.risk for loss in dd.risk.engineering.loss)
 
-    #= ====== =#
+#= ====== =#
 #  limits  #
 #= ====== =#
 dyexp["limits.model[:].cleared"] =

From 4474394288118630a8a405f3011b298f84fe01a0 Mon Sep 17 00:00:00 2001
From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com>
Date: Wed, 5 Mar 2025 14:14:21 -0700
Subject: [PATCH 6/6] Extract and objective functions for risk

---
 src/extract/extract.jl    | 3 +++
 src/extract/objectives.jl | 3 +++
 2 files changed, 6 insertions(+)

diff --git a/src/extract/extract.jl b/src/extract/extract.jl
index 6539aab6..7b728f73 100644
--- a/src/extract/extract.jl
+++ b/src/extract/extract.jl
@@ -380,6 +380,9 @@ function update_ExtractFunctionsLibrary!()
     ExtractLibFunction(:costing, :blanket_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"blanket") / dd.costing.cost_direct_capital.cost)
     ExtractLibFunction(:costing, :cryostat_of_total, "%", dd -> 100 * select_direct_captial_cost(dd,"cryostat") / dd.costing.cost_direct_capital.cost)
 
+    ExtractLibFunction(:risk, :total_engineering_risk, "\$M", dd -> dd.risk.engineering.risk)
+    ExtractLibFunction(:risk, :total_plasma_risk, "\$/kWh", dd -> dd.risk.plasma.risk)
+ 
     #! format: on
 
     return EFL
diff --git a/src/extract/objectives.jl b/src/extract/objectives.jl
index d43b0b0f..7fa28f92 100644
--- a/src/extract/objectives.jl
+++ b/src/extract/objectives.jl
@@ -97,6 +97,9 @@ function update_ObjectiveFunctionsLibrary!()
     ObjectiveFunction(:min_βn, "", dd -> dd.equilibrium.time_slice[].global_quantities.beta_normal, -Inf)
     ObjectiveFunction(:min_R0, "m", dd -> dd.equilibrium.time_slice[].boundary.geometric_axis.r, -Inf)
     ObjectiveFunction(:max_zeff, "", dd -> @ddtime(dd.summary.volume_average.zeff.value), Inf)
+    ObjectiveFunction(:min_engineering_risk, "\$M", dd -> dd.risk.engineering.risk, -Inf)
+    ObjectiveFunction(:min_plasma_risk, "\$/kWh", dd -> dd.risk.plasma.risk, -Inf)
+
     #! format: on
     return ObjectiveFunctionsLibrary
 end