Use consistent form in constraint equations

[jonmaddock]

Following an investigation into upper limit constraint behaviour, I'm proposing to change the form of the upper limit inequality constraints from:
$$c_{max} = 1 - \frac{x_{max}}{x}$$
to
$$c_{max} = \frac{x}{x_{max}} - 1$$
This means that lower and upper limit constraints are violated proportionately (e.g. linearly), and the constraint values for upper limit constraints are now much simpler to interpret. It also removes a discontinuity in upper limit constraints when the constrained parameter ($x$) passes through 0.

In particular, I'd be grateful for reviewers to carefully check my changes to constraints 12, 15 and 86, which I believe were incorrectly (or at least inconsistently) defined. I would also like 76 to be checked, as I'm not sure I've got this the right way around.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 31.66%. Comparing base (2e441df) to head (10c4b3f).

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #3565   +/-   ##
=======================================
  Coverage   31.66%   31.66%           
=======================================
  Files          86       86           
  Lines       20198    20198           
=======================================
  Hits         6396     6396           
  Misses      13802    13802           

☔ 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.

[mkovari]

76 looks OK to me:

OLD:   tmp_cc = 1.0D0 - fnesep * nesep_crit/nesep
NEW:   tmp_cc = nesep / nesep_crit - 1.0D0 * fnesep

In both cases nesep_crit is an upper limit for nesep, if fnesep is less than or equal to 1.

[mkovari]

86 also looks OK to me:

OLD      tmp_cc = 1.0D0 - t_turn_tf / ( f_t_turn_tf * t_turn_tf_max )
NEW      tmp_cc = t_turn_tf / t_turn_tf_max - 1.0D0 * f_t_turn_tf

In both cases t_turn_tf_max is an upper limit for t_turn_tf , if f_t_turn_tf is less than or equal to 1.

[jonmaddock]

5% regression test differences

I've tried to help understand the regression test differences as a result of these changes. 6 failed, 1 passed.

large_tokamak_no_f

297 differences. Inequality values dominate unsurprisingly, but there are some significant other differences; some of the largest are shown below. A different solution is being found (as in our previous discussion @mkovari ): beta has increased by nearly 20%.

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
pcurrentdrivemw     	      7.53	      37.8	    402.03
auxiliary_cd        	  3.46e+05	  1.67e+06	    381.96
aux_current_fraction	    0.0205	    0.0936	    356.46
divlife_cal         	      3.52	      12.2	    246.50
divlife             	       4.4	      15.2	    246.50
boundl(103)         	      93.5	    0.0908	    -99.90
pfc1t4              	 -1.15e+06	 -6.33e+05	     45.06
pfc0t4              	  1.26e+06	  1.72e+06	     37.14
p_plasma_outer_rad_mw	       141	       193	     36.85
c2231               	       248	       338	     36.69
c223                	       248	       338	     36.69
echwpow             	       165	       226	     36.69
pheatingmw          	       165	       226	     36.69
pinjmw.             	      82.5	       113	     36.69
pinjemw             	      82.5	       113	     36.69
echpwr              	      82.5	       113	     36.69
pinjht              	      82.5	       113	     36.69
pinjmw              	      82.5	       113	     36.69
pinjwp              	       165	       226	     36.69
nd_impurities       	  4.85e+16	   6.6e+16	     36.00
p_plasma_rad_mw     	       236	       314	     33.16
pflux_fw_rad_max_mw 	     0.616	      0.82	     33.16
pflux_fw_rad_mw     	     0.185	     0.246	     33.16
fimp(13)            	  0.000588	   0.00077	     30.82
xcm020              	      1.55	      2.03	     30.82
itvar020            	  0.000588	   0.00077	     30.82
pfc4t1              	  4.07e+05	   5.2e+05	     27.95
pfc5t1              	  4.07e+05	   5.2e+05	     27.95
pcoreradmw.         	      94.7	       121	     27.67
p_plasma_inner_rad_mw	      94.7	       121	     27.67
radiation_power_subtracted_from_plasma_power_balance_(mw)	      94.7	       121	     27.67
whtconsc            	  3.16e+03	  2.43e+03	    -23.20
beta_thermal_toroidal	    0.0341	    0.0411	     20.52
beta_mcdonald       	    0.0389	    0.0468	     20.37
beta_toroidal       	    0.0389	    0.0468	     20.37
beta_thermal        	    0.0331	    0.0396	     19.70
beta                	    0.0377	    0.0451	     19.55
itvar004            	    0.0377	    0.0451	     19.55
...

spherical_tokamak_once_through

Only 1 difference >5%, but a dramatic one. In this case , eq_con046 is "Ip/Irod_upper_limit_normalised_residue": it is in fact an upper limit, which has had its form changed. This would result in a different solution (were it an optimisation problem), as this violation is now above the solver tolerance, but it should be compensated for by its f-value adjusting accordingly.

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
eq_con046           	 -2.33e-15	     -1.29	-55124764047814096.00

helias_5b

Again, only 1 difference. ineq_con024 = "Beta_upper_limit_normalised_residue": it's an upper limit, so we might expect it to change, or a slightly different solution is being found.

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
ineq_con024         	     0.101	    0.0915	     -9.15

large_tokamak_once_through

12 diffs, all of them inequality constraints, as expected.

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
ineq_con015         	     0.433	     0.764	     76.37
ineq_con030         	      1.64	     0.621	    -62.11
ineq_con024         	      0.21	    0.0869	    -58.69
ineq_con065         	      1.27	      0.56	    -55.99
ineq_con008         	     0.961	      0.49	    -49.01
ineq_con027         	     0.149	    0.0779	    -47.79
ineq_con009         	      0.84	     0.457	    -45.66
ineq_con026         	    0.0554	    0.0315	    -43.15
ineq_con033         	   0.00155	   0.00101	    -35.10
ineq_con005         	   0.00139	   0.00167	     19.83
ineq_con032         	     0.175	     0.149	    -14.88
ineq_con025         	     0.144	     0.126	    -12.61

stellarator_helias_once_through

7 diffs, again all inequalities.

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
ineq_con017         	      5.41	     0.844	    -84.41
ineq_con065         	      5.03	     0.834	    -83.41
ineq_con018         	      3.18	     0.609	    -80.86
ineq_con032         	      2.53	     0.717	    -71.67
ineq_con067         	     0.431	     0.301	    -30.13
ineq_con034         	     0.195	     0.163	    -16.30
ineq_con024         	     0.153	     0.133	    -13.25

large_tokamak

62 diffs. This is an equality-only problem, so some equality diffs are expected, but a different solution is being found too (e.g. auxiliary_cd).

Variable            	       Ref	       New	  % Change
------------------------------------------------------------
eq_con015           	 -1.44e-09	  5.73e-07	  39971.46
eq_con036           	 -7.79e-07	 -1.52e-05	  -1853.83
eq_con026           	   1.6e-08	 -2.07e-07	  -1397.56
eq_con002           	   1.9e-08	   2.8e-07	   1376.24
eq_con033           	 -2.92e-07	 -4.28e-06	  -1364.68
eq_con013           	 -2.48e-07	  2.33e-06	   1041.07
eq_con030           	 -6.68e-08	  5.93e-07	    986.76
eq_con062           	  5.99e-08	 -5.18e-07	   -965.45
eq_con060           	  6.43e-08	 -4.96e-07	   -870.71
eq_con001           	  1.37e-07	  1.17e-06	    748.16
eq_con068           	  1.35e-07	   9.6e-07	    612.07
eq_con027           	  1.84e-08	 -8.01e-08	   -536.34
eq_con072           	 -6.83e-08	 -3.04e-07	   -345.78
eq_con008           	  1.11e-07	 -1.94e-07	   -275.10
eq_con016           	 -4.21e-07	  5.68e-07	    234.80
eq_con025           	 -9.74e-09	 -2.87e-08	   -194.40
eq_con035           	 -2.72e-08	  9.33e-09	    134.33
pcurrentdrivemw     	      4.18	      9.57	    128.98
auxiliary_cd        	  1.94e+05	  4.43e+05	    128.03
aux_current_fraction	    0.0117	    0.0265	    127.09
eq_con024           	  5.14e-08	 -1.04e-08	   -120.14
eq_con011           	  7.91e-10	 -1.41e-10	   -117.86
n_cycle             	         0	       256	    100.00
boundl(103)         	      95.2	    0.0944	    -99.90
itvar034            	      1.71	     0.582	    -66.04
xcm034              	      1.71	     0.582	    -66.04
...

Summary

Due to having equality-only, mixed equality-inequality and model-evaluation (non-optimising) regression tests, we can see the impact of these changes on optimised solutions as well as model evaluations. Whilst it isn't surprising that the inequality constraint values change in the model evaluation case, it is interesting to see that significantly different optimised solutions are being found in some cases as a result of changing just the constraint values, not their limits; the value of the constraint functions affects the solution, not just the position of their limits. As @mkovari suggested previously, this could be because these problems are underdetermined; as a result of these changes the solver is now finding another non-unique solution instead. We should probably pay more attention to the possibility of non-unique solutions in future.

[jonmaddock]

@ajpearcey @mkovari @stuartmuldrew can I have a review for this please?

[ajpearcey]

Looks to be ready to merge to me.

[mkovari]

Sorry, @jonmaddock , I managed to delete my comment relating to the sign of the central solenoid current. I think once this is checked the changes can be merged, but it would be good in the near future to see the results once the lower bound for rmajor is lowered.

This was my previous comment:
The correct sign here depends on the sign of vstot and vs_plasma_total_required.

[jonmaddock]

Before these changes, evaluating `large_tokamak_once_through_IN.DAT:

 vs_cs_pf_total_pulse =  -577.22422897684089     
 vs_plasma_total_required =   557.98626457279624     
 con12_tmp_cc =  -3.4477487396170181E-002

After these changes:

 vs_cs_pf_total_pulse =  -577.22422897684089     
 vs_plasma_total_required =   557.98626457279624     
 con12_tmp_cc =   2.0344774873961704  

The constraint has changed sign, which is incorrect. My modification to this was wrong, so thank you @mkovari for spotting this! Correcting for the negative sign of vs_cs_pf_total_pulse, this produces:

 vs_cs_pf_total_pulse =  -577.22422897684089     
 vs_plasma_total_required =   557.98626457279624     
 tmp_cc =  -3.4477487396170181E-002

... the same as originally. I've pushed this correction. If you're happy with this @mkovari , please can you approve the review and I'll merge. Thanks for spotting this.

[timothy-nunn]

Can you double check this constraint.

I have run a once-through at the initial ST design point on main and on this branch. The only constraint with a sign change is 046:

On main its feasible at the initial point:

Ip/Irod_upper_limit_normalised_residue___________________________________ (eq_con046)____________________ 9.20229325178258950e-02

but on this branch it is not:

Ip/Irod_upper_limit_normalised_residue___________________________________ (eq_con046)____________________ -1.30279815951349276e+00

Could this be because cratmx is on the top of the constraint when I thought it should be on the bottom? Maybe this is why the ST regression test is failing?

[jonmaddock]

You are absolutely right about this Tim, well spotted! I've changed this. However, when I run this on main, I get:

 Ip/Irod_upper_limit_normalised_residue___________________________________ (eq_con046)____________________ -1.88737914186276612e-15

i.e. slightly violated. After my (corrected) change I get:

 Ip/Irod_upper_limit_normalised_residue___________________________________ (eq_con046)____________________ -9.99200722162640886e-16

again violated by a tiny amount. The version of main I ran on is the same that I rebased this branch on.

However, the ST regression test now solves again, brilliant!

[timothy-nunn]

Glad to hear it, as long as main and this branch have the same sign that is all that matters!

[timothy-nunn]

Getting this warning too. Please check this is not causing a problem.

  889 |       tmp_cc =  1.0D0 - fvs * -vs_cs_pf_total_pulse/vs_plasma_total_required
      |                               1
Warning: Extension: Unary operator following arithmetic operator (use parentheses) at (1)

[jonmaddock]

Quite right, that was careless! I've enclosed it.

[ajpearcey]

I should actually formally approve this, following my review.