Skip to content

Conversation

@timothy-nunn
Copy link
Collaborator

@timothy-nunn timothy-nunn commented Oct 22, 2025

Removes F-values from PROCESS

  • Add a test for the constraints that runs and checks no attribute/type errors occur (fixed a few bugs that arose before the removal)
  • Removed f-values from constraints, iteration variables, inputs, scan variables, and their mention in the docs. The follow f-values have been kept but have been removed from their ability to iterate:
    • fl_h_threshold
    • fiooic
    • fjohc
    • fjohc0
    • fdene
    • fradpwr
  • Removed f-values from all regression tests
    • st_regression.IN.DAT used some weird f-value bounds so I have updated the constraints (see comment below).
  • Added an additional constraint (22) to enforce L-mode (constraint 15 remains to enforce H-mode).

@timothy-nunn timothy-nunn linked an issue Oct 22, 2025 that may be closed by this pull request
@timothy-nunn timothy-nunn self-assigned this Oct 30, 2025
@codecov-commenter
Copy link

codecov-commenter commented Nov 5, 2025

Codecov Report

❌ Patch coverage is 80.85106% with 9 lines in your changes missing coverage. Please review.
✅ Project coverage is 44.20%. Comparing base (8a4d67c) to head (f9a5a9b).
⚠️ Report is 3 commits behind head on main.

Files with missing lines Patch % Lines
process/pfcoil.py 0.00% 4 Missing ⚠️
process/constraints.py 80.00% 3 Missing ⚠️
process/superconducting_tf_coil.py 0.00% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #3947      +/-   ##
==========================================
- Coverage   46.30%   44.20%   -2.11%     
==========================================
  Files         123      123              
  Lines       28961    30581    +1620     
==========================================
+ Hits        13411    13518     +107     
- Misses      15550    17063    +1513     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@timothy-nunn timothy-nunn force-pushed the 3936-remove-f-values branch 3 times, most recently from 1bd708f to 6302f99 Compare November 6, 2025 13:45
Comment on lines 2685 to 2692
ixc = 46
fp_hcd_injected_max = 1.0
boundl(46) = 0.6
boundu(46) = 1.5
* DESCRIPTION: f-value for injected power variation range
* JUSTIFICATION: Setup to allow the injected power to vary

p_hcd_injected_max = 150.0
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When fp_hcd_injected_max = 1.5: $\frac{225.0}{150.0}-1.5 = 0$ therefore the new bound should be $225.0$

*--------------------------Numerics Variables---------------------------*
*‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾*

epsvmc = 1.0E-11
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was much lower than other regression tests

@timothy-nunn timothy-nunn marked this pull request as ready for review November 13, 2025 16:54
@timothy-nunn
Copy link
Collaborator Author

How to update an input file to not use f-values

  1. Ensure all equality constraint definitions (icc = ...) occur before any inequality constraints.
  2. Set the variable neqns equal to the number of equality constraints.
  3. Remove obsolete f-values from the input file.
  4. Remove iteration variables that corresponded to f-values.

In 99% of cases, this will allow your file to run without f-values and without treating inequality constraints as equality. However, if the input file specified weird bounds on the f-value iteration variables, you will need to modify the constraint bound (max/min) to account for this change.

Copy link
Collaborator

@mkovari mkovari left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See comments



!!! Info "Constraints"
A full list of constraints is given on the variable description page in the row labelled
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need to update this in due course.

OBS_VARS_HELP.update(
dict.fromkeys(
fvalues_list,
"F-values are no longer supported, please use inequality constraints",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When the inequality constraints are documented we could add a link to the docs here.

numerics.boundu[37]
* pfcoil_variables.j_cs_critical_flat_top_end
# numerics.boundu[fjohc] * (which should have been 1)
pfcoil_variables.j_cs_critical_flat_top_end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this can ever be false:
abs(pfcoil_variables.j_cs_flat_top_end) > 0.99e0 * abs(pfcoil_variables.j_cs_critical_flat_top_end)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Keep this warning (cslimit).
Keep fjohc and fjohc0 as variables (not iteration) and include them in the constraint.
Use them instead of the bounds in this line of code.

"CS not using max current density: further optimisation may be possible"
)

# Check whether CS coil currents are feasible from engineering POV
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally these warnings should be recreated but it's no big deal.


# The operation current density weighted with the global iop/icrit fraction
lhs[:] = constraint_variables.fiooic * jcrit_vector
lhs[:] = jcrit_vector
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't really correct anymore, I suspect, and the same may be true for the code below.

Copy link
Collaborator Author

@timothy-nunn timothy-nunn Nov 25, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the recourse here then? Should fiooic be added back? Or is there a way to calculate it from other variables in PROCESS (maybe the mentioned iop/icrit)?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I really don't know anything about the stellarator code. I don't know what lhs and rhs are.


This constraint can be activated by stating `icc = 15` in the input file.

The scaling value `fl_h_threshold (ixc=103)` can be varied to set the required margin around the threshold.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fl_h_threshold should be kept and renamed. This new way of using it should be documented here.

Copy link
Collaborator Author

@timothy-nunn timothy-nunn Nov 25, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have updated the description and docs to document how I think fl_h_threshold is working now. We should probably go through this at some point (with some example input files) to ensure its working as before/expected.

In particular, we should carefully check how it was used to enforce L-mode/H-mode before.

@timothy-nunn
Copy link
Collaborator Author

timothy-nunn commented Nov 25, 2025

Hi @ajpearcey @mkovari

Thank you both for the reviews. I have resolved everything that I have actioned. The following comments are still outstanding:

@timothy-nunn
Copy link
Collaborator Author

timothy-nunn commented Nov 25, 2025

@chris-ashe has pointed me to f_c_tf_turn_operating_critical which is possibly the same as fiooic.

Looking at an MFile, they have very similar values:

fiooic___________________________________________________________________ (itvar024)_____________________ 7.28529047768881166e-01 
Actual_current_/_critical_current________________________________________ (f_c_tf_turn_operating_critical)_ 7.28529093232043401e-01 OP 

So I can replace the warning and its use in stellarator with f_c_tf_turn_operating_critical.

This will not work for the Stellarator (#3947 (comment)) because f_c_tf_turn_operating_critical is not set, nor are a lot of variables related to the TF critical current (@ym1906 any ideas since you're probably most familiar with the Stellarator models? I could add an input in that is specifically a margin for the Stellarator model that would not iterate and would default to say 0.7).

@mkovari
Copy link
Collaborator

mkovari commented Nov 26, 2025

There are two separate constraint equations for enforcing the L-H threshold: icc = 15 and 73.
In both cases it would be tidier to have separate constraints for L and H mode.
In any case there needs to be a margin defined by a parameter (not an iteration variable) to replace fl_h_threshold.

@timothy-nunn
Copy link
Collaborator Author

There are two separate constraint equations for enforcing the L-H threshold: icc = 15 and 73. In both cases it would be tidier to have separate constraints for L and H mode. In any case there needs to be a margin defined by a parameter (not an iteration variable) to replace fl_h_threshold.

@ajpearcey @chris-ashe @jonmaddock do we have any thoughts on this? I'm not sure we can use fl_h_threshold as set out here.
image

I think this requires the margin (f-value) to be iterating.


# The operation current density weighted with the global iop/icrit fraction
lhs[:] = constraint_variables.fiooic * jcrit_vector
lhs[:] = jcrit_vector
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I really don't know anything about the stellarator code. I don't know what lhs and rhs are.

numerics.boundu[37]
* pfcoil_variables.j_cs_critical_flat_top_end
# numerics.boundu[fjohc] * (which should have been 1)
pfcoil_variables.j_cs_critical_flat_top_end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Keep this warning (cslimit).
Keep fjohc and fjohc0 as variables (not iteration) and include them in the constraint.
Use them instead of the bounds in this line of code.

"j_cs_pulse_start / j_cs_critical_pulse_start shouldn't be above 0.7 "
"for engineering reliability"
)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These warnings are fine. Remember the constraints on current density may not be applied (if T margin is used instead).

"f_c_tf_turn_operating_critical shouldn't be above 0.7 for engineering reliability"
)

po.ovarre(
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This warning is fine.

args : output structure : residual error; constraint value;
fiooic: f-value for TF coil operating current / critical
j_tf_wp_critical: critical current density for winding pack (A/m2)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We still need fiooic:
The current density margin can be set using fiooic.

This constraint can be activated by stating `icc = 5` in the input file.

The value of `i_density_limit` can be set to apply the relevant limit . The scaling value `fdene` can be varied also.
The value of `i_density_limit` can be set to apply the relevant limit. The variable `fdene` can be set to scale the constraint bound: `nd_plasma_electrons_vol_avg` / `nd_plasma_electrons_max` <= `fdene`.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i_density_limit isnt always limits that apply to the volume averaged density.

$$

**Therefore it is recommended to always use `icc = 15` if trying to simulate a plasma scenario specifically in L or H-mode**
Again, `fl_h_threshold` can be used to add a margin to the constraint. For example, `fl_h_threshold = 1.25` ensures that `p_plasma_separatrix_mw` is at least 25% less than the threshold power.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should it not be 25% more?

******************
nsweep = 17 * b_tf_inboard_max, maximum peak toroidal field
isweep = 6
sweep = 10.5, 10.4, 10.3, 10.2, 10.1, 10.0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we know how close the solution is initially to this inboard leg limit? Would it be better to have a scan example that directly changes a variable instead of a limit?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think for this PR I will leave it but we should probably create an issue to make a more meaningful example

)


@ConstraintManager.register_constraint(22, "MW", ">=")
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the sake of numerical sanity should this not be further down?

@timothy-nunn
Copy link
Collaborator Author

Chris also recommends having two separate margins for constraint 15 and 22

@timothy-nunn
Copy link
Collaborator Author

timothy-nunn commented Dec 15, 2025

Hi All, I have added in the two separate variables and updated the docstrings/docs accordingly. Please check my work carefully to ensure this new way of constraining on the LH threshold makes sense.

I have had to make some changes in physics.py because fl_h_threshold was used there, I am not sure if what I have done is correct?


$$
1.0 - \mathtt{fl\_h\_threshold} \times \frac{\overbrace{\mathtt{p\_l\_h\_threshold\_mw}}^{\text{Power from scaling}}}{\mathtt{p_plasma_separatrix_mw}}
\mathtt{p\_plasma\_separatrix\_mw} \ge \mathtt{h\_mode\_threshold\_margin} \times \underbrace{\mathtt{p\_l\_h\_threshold\_mw}}_{\text{Power from scaling}}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For writing variable names as equations I learned its always better to use \texttt{} as it prevents having to use \_ for the spacing. Also means the rendering doesn't break when the name is changed.

@CoronelBuendia
Copy link
Contributor

Hi @timothy-nunn, hi all,

Apologies in advance for what is probably not a particularly actionable comment, and which may lead only to a long and inconclusive dive into the literature… but after this, I promise to forever hold my peace!

The existence of f-values in PROCESS has intrigued me for a long time. I’ve never been entirely confident in linear programming or convex optimisation theory, so I’m still a little hesitant to bring this up and may not be using the right terminology.

At first glance, converting all inequality constraints into equalities seems completely counterintuitive — it makes the problem harder, not easier. But there are two thoughts I can’t shake:

  • The original PROCESS authors (whose code was deeply intertwined with VMCON) really knew their maths and optimisation — they must have done this for a reason.
  • PROCESS’s optimisation mode has, on average and across versions, been surprisingly robust and reliable.

The f-value formulation reminds me of duality in optimisation: the idea that a linear program can be reformulated in another equivalent form. This “standard reformulation” can sometimes be easier to solve, counterintuitive as that sounds, because gradient-based optimisers often handle linear equalities more gracefully — they can be treated similarly to bounds. My hunch is that the old VMCON was actually designed to work best with this kind of formulation. It’s all KKT and Lagrange multiplier territory — though I’m not sure how the new VMCON handles this compared to the old one.

Perhaps the simplest way to put my doubts to rest would be to show that an EU-DEMO (non-LTS) regression test still converges in the same way.

The other potential uses of f-values I see are:

  • as margins (heavily discussed here and I have nothing to add)
  • as fudge factors when trying to go back and forth between PROCESS and something at higher fidelity.

The first case is already well covered in this PR, and the second could be implemented more straightforwardly on the user side — simply by adjusting constraint values. So, in my view, neither case really justifies the existence of f-values as a built-in feature.

@mkovari
Copy link
Collaborator

mkovari commented Jan 14, 2026

Thanks for continuing to educate us...

I like this line in the Wikipedia article:

According to George Dantzig, the duality theorem for linear optimization was conjectured by John von Neumann immediately after Dantzig presented the linear programming problem.

  • First, the slack variables (f-values) are not the Lagrange multipliers.
  • (b) VMCON was fully equipped with the facility to handle equality and inequality constraints, and I trust Powell, Crane et al more than the people who wrote PROCESS.
  • iii. I seem to remember that the justification for the use of slack variables in PROCESS was to make it possible to use the non-optimising "HYBRID" solver, which didn't allow constraints at all.
  • Finally, PROCESS’s optimisation mode has never been robust or reliable.
  • According to Wikipedia, Slack variables are used in particular in linear programming. We are not doing linear programming.

\texttt{p_plasma_separatrix_mw} \ge \texttt{f_h_mode_margin} \times \underbrace{\texttt{p_l_h_threshold_mw}}_{\text{Power from scaling}}
$$

For example, `f_h_mode_margin = 1.2` ensures that `p_plasma_separatrix_mw` is at least $1.2\times$ greater than the threshold power `p_l_h_threshold_mw`.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Has the constraint been updated? to stop discontinuities in the constraint fl_h_threshold is the inverse. For example to have p_plasma_separatrix_mw above 1.2 times p_l_h_threshold_mw, you need f_h_mode_margin = 0.833

Copy link
Collaborator

@ajpearcey ajpearcey Jan 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

okay ignore me, looks like you've already done this.

@timothy-nunn
Copy link
Collaborator Author

Hi @CoronelBuendia,

Thank you for weighing in on this, you raise important points and share a lot of the concerns I have had: why were f-values ever added (among other things).

I have gone digging back on our GitLab and have found the following written in a file process.tex that existed at least 11 years ago, and probably before:

Limit equations are not restricted to optimisation mode. In non-optimisation
mode, the iteration variables are not bounded, but the f-values can still be
used to provide information about how calculated values compare with limiting
values, without having to change the characteristics of the device being
benchmarked to find a solution.

It is for this reason that all the constraint equations used in \process\ are
formulated as equalities, despite the fact that equation solver \vmcon\ can
solve for inequalities as well. The use of f-values precludes this need, and
allows the non-optimising equation solver HYBRD to use the same constraint
equations.

I therefore agree with @mkovari that f-values were only ever added for the HYBRD solver that could not support inequality constraints. We no longer use HYBRD and f-values are now only (mis)used to apply margin to constraints. In most cases (as you say), this can be done by simply modifying the constraint bound which is much more transparent; in some cases we have kept margin variables but not allowed them to iterate anymore.

We have conducted several convergence studies over the last year or two and there is little difference in the convergence of large_tokamak and large_tokamak_nof when randomly varying the initial iteration variables by +/- 50%. This tells me that the removal of f-values neither helps nor hinders the solver (which is surprising because we roughly halve the dimensionality of the problem!).

@CoronelBuendia
Copy link
Contributor

@mkovari Ha! Educating is a generous word to use non-sarcastically for that fuzzy drivel I spewed... alternatively - good sarcasm! As I said, I'm shaky on the foundations. LP and QP are very similar though... only the form of the objective function differs. Note that minimising R_0 is technically an LP, as $R_0 \in \mathbf{x}$, even if under the hood VMCON is doing sequential QPs.

I think that, actually, for the complexity of the optimisation problems typically solved in PROCESS, its optimisation mode is remarkably reliable. Which is not to say that it works every time or even most of the time; it's just often a non-trivial problem with models that fight each other (sometimes in ways that fail) but that isn't really the fault of the optimisation in itself. Perhaps I just have too little faith in it to begin with.. so I've often been pleasantly surprised.

@mkovari @timothy-nunn I was not aware that the primary use case of f-values was for HYBRD. This clarifies a lot. I'm glad this conversation will live in the history (on Github at least). I don't really want to create work, but this might be a good one for a short summary ADR, so that it lives in the code, too.

All that said, with regards to the regression test, I think that large_tokamak is a really bad yardstick. The formulation is strange, the options are weird, important and historically problematic constraints are not used, and the optimal result is a lower bound of an iteration variable, meaning that some constraints are irrelevant. The hardest gradient-descent convergence paths are when multiple non-linear constraints are pulling in different directions in a way that affects the objective function value. Although, in fairness, this does depend on the convergence criterion[ia] used.

And now, I will forever hold my peace :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Remove f-values

7 participants