Optimum solutions comparison notebook¶

A Jupyter notebook to demonstrate how changing input parameters changes the optimum solution found by PROCESS.

NOTE We run the examples in a temporary directory so all the inputs are copied there and the outputs contained there before the directory is removed when the example has finished running. This keeps the examples directory tidy and does not permanently modify any data files. The use of temporary directories is not needed for regular use of PROCESS.

This notebook demonstrates how the optimum solution found by PROCESS changes as we vary input parameters. We will use the large tokamak example input file to do this. The figure of merit for this example is to minimise the major radius, rmajor.

We use the functionality from plot_solutions.py and from plot_proc.py to demonstrate this.

These tools plot the solution vectors (i.e. final values of optimisation parameters) for different runs of PROCESS. This allows visual comparisons of different solution points.

It can use different intra-solution optimisation parameter normalisations (e.g. initial value, parameter range) and inter-solution normalisations (e.g. normalise to a certain solution).

Known Limitations¶

• The solution vectors (optimisation parameter values at the solution) currently plotted are normalised to the initial point (from the IN.DAT) of each solution: each element of the vector is the $x_{final}/x_{initial}$, the xcmxxx values in the MFILE.DAT. This allows all optimisation parameters to be plotted on the same axis, showing the relative changes from their initial values across multiple solutions.
• Solutions being plotted together must also have the same optimisation parameters.
• The solutions plotted in this example are fictitious.

Setup¶

First we need to generate the necessary MFILES to be used in the rest of this notebook.

In [1]:

%load_ext autoreload
%autoreload 2

import shutil
import tempfile
from pathlib import Path

from process.core.io.plot_solutions import (
    RunMetadata,
    plot_mfile_solutions,
)
from process.core.repository import get_process_root
from process.main import SingleRun

working_dir = Path.cwd()

# Define input file name relative to project dir, then copy to temp dir
script_dir = Path("__file__").parent.resolve()

data_dir = get_process_root() / "../examples/data/"
input_file = data_dir / "large_tokamak_IN.DAT"
input_file2 = data_dir / "large_tokamak_varied_min_net_electric_IN.DAT"

%load_ext autoreload %autoreload 2 import shutil import tempfile from pathlib import Path from process.core.io.plot_solutions import ( RunMetadata, plot_mfile_solutions, ) from process.core.repository import get_process_root from process.main import SingleRun working_dir = Path.cwd() # Define input file name relative to project dir, then copy to temp dir script_dir = Path("__file__").parent.resolve() data_dir = get_process_root() / "../examples/data/" input_file = data_dir / "large_tokamak_IN.DAT" input_file2 = data_dir / "large_tokamak_varied_min_net_electric_IN.DAT"

In [2]:

# Copy the file to avoid polluting the project directory with example files
temp_dir = tempfile.TemporaryDirectory()
working_dir = Path(temp_dir.name)
input_path = Path(temp_dir.name) / "large_tokamak_IN.DAT"
input_path2 = Path(temp_dir.name) / "large_tokamak_varied_min_net_electric_IN.DAT"
shutil.copy(input_file, input_path)
shutil.copy(input_file2, input_path2)

# Copy the file to avoid polluting the project directory with example files temp_dir = tempfile.TemporaryDirectory() working_dir = Path(temp_dir.name) input_path = Path(temp_dir.name) / "large_tokamak_IN.DAT" input_path2 = Path(temp_dir.name) / "large_tokamak_varied_min_net_electric_IN.DAT" shutil.copy(input_file, input_path) shutil.copy(input_file2, input_path2)

Out[2]:

PosixPath('/tmp/tmpww7a3bhb/large_tokamak_varied_min_net_electric_IN.DAT')

In [3]:

# Run process on these input files in a temporary directory
single_run = SingleRun(input_path.as_posix())
single_run.run()
single_run = SingleRun(input_path2.as_posix())
single_run.run()

# Run process on these input files in a temporary directory single_run = SingleRun(input_path.as_posix()) single_run.run() single_run = SingleRun(input_path2.as_posix()) single_run.run()

The IN.DAT file does not contain any obsolete variables.
 
**************************************************************************************************************
************************************************** PROCESS ***************************************************
************************************** Power Reactor Optimisation Code ***************************************
**************************************************************************************************************
 
Version : 3.3.1.dev59+g87727f966
Git Tag : v3.3.0-59-g87727f96
Git Branch : main
Date : 25/03/2026 UTC
Time : 16:04
User : runner
Computer : runnervm46oaq
Directory : /home/runner/work/PROCESS/PROCESS
Input : /tmp/tmpww7a3bhb/large_tokamak_IN.DAT
Run title : Generic large tokamak
Run type : Reactor concept design: Pulsed tokamak model, (c) UK Atomic Energy Authority
 
**************************************************************************************************************
 
Equality constraints : 3
Inequality constraints : 23
Total constraints : 26
Iteration variables : 19
Max iterations : 200
Figure of merit : +1  -- minimise major radius          
Convergence parameter : 1e-07
 
**************************************************************************************************************

/home/runner/work/PROCESS/PROCESS/process/core/init.py:92: UserWarning: Lower limit of volume averaged electron temperature (temp_plasma_electron_vol_avg_kev) has been raised to ensure temp_plasma_electron_vol_avg_kev > temp_plasma_pedestal_kev
  check_process(inputs)
/home/runner/work/PROCESS/PROCESS/process/core/init.py:92: UserWarning: temp_cs_superconductor_margin_min and tmargmin should not both be specified in IN.DAT temp_cs_superconductor_margin_min has been ignored
  check_process(inputs)

/home/runner/work/PROCESS/PROCESS/process/models/physics/bootstrap_current.py:1194: RuntimeWarning: divide by zero encountered in scalar divide
  * (nd_plasma_pedestal_electron / n_greenwald) ** -0.174

1 | Convergence Parameter: 7.936E-01

/home/runner/work/PROCESS/PROCESS/process/models/physics/current_drive.py:2154: RuntimeWarning: invalid value encountered in scalar divide
  * (c_hcd_driven / p_hcd_injected)

2 | Convergence Parameter: 1.831E-01

3 | Convergence Parameter: 2.837E-02

4 | Convergence Parameter: 7.247E-02

5 | Convergence Parameter: 2.120E-02

6 | Convergence Parameter: 8.867E-03

7 | Convergence Parameter: 1.631E-02

8 | Convergence Parameter: 2.932E-03

9 | Convergence Parameter: 2.889E-03

10 | Convergence Parameter: 6.735E-05

11 | Convergence Parameter: 4.008E-04

12 | Convergence Parameter: 3.389E-03

13 | Convergence Parameter: 2.045E-02

14 | Convergence Parameter: 1.039E-03

15 | Convergence Parameter: 1.224E-07

16 | Convergence Parameter: 6.642E-10

17 | Convergence Parameter: 2.462E-09

 
************************************* PROCESS found a feasible solution **************************************
 

******************************************** Errors and Warnings *********************************************
(/home/runner/work/PROCESS/PROCESS/process/models/physics/bootstrap_current.py:1398) Diamagnetic fraction is more than 1%, but not calculated. Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1

(/home/runner/work/PROCESS/PROCESS/process/models/tfcoil/base.py:255) dr_tf_plasma_case too small to accommodate the WP, forced to minimum value
 
******************************************* End of PROCESS Output ********************************************
 
The IN.DAT file does not contain any obsolete variables.
 
**************************************************************************************************************
************************************************** PROCESS ***************************************************
************************************** Power Reactor Optimisation Code ***************************************
**************************************************************************************************************
 
Version : 3.3.1.dev59+g87727f966
Git Tag : v3.3.0-59-g87727f96
Git Branch : main
Date : 25/03/2026 UTC
Time : 16:05
User : runner
Computer : runnervm46oaq
Directory : /home/runner/work/PROCESS/PROCESS
Input : /tmp/tmpww7a3bhb/large_tokamak_varied_min_net_electric_IN.DAT
Run title : Generic large tokamak
Run type : Reactor concept design: Pulsed tokamak model, (c) UK Atomic Energy Authority
 
**************************************************************************************************************
 
Equality constraints : 3
Inequality constraints : 23
Total constraints : 26
Iteration variables : 19
Max iterations : 200
Figure of merit : +1  -- minimise major radius          
Convergence parameter : 1e-07
 
**************************************************************************************************************

1 | Convergence Parameter: 2.627E-01

2 | Convergence Parameter: 5.821E-02

3 | Convergence Parameter: 3.613E-02

4 | Convergence Parameter: 1.024E-03

5 | Convergence Parameter: 7.297E-10

6 | Convergence Parameter: 2.094E-09

 
************************************* PROCESS found a feasible solution **************************************
 

******************************************** Errors and Warnings *********************************************
(/home/runner/work/PROCESS/PROCESS/process/models/physics/bootstrap_current.py:1398) Diamagnetic fraction is more than 1%, but not calculated. Consider using i_diamagnetic_current=2 and i_pfirsch_schluter_current=1

(/home/runner/work/PROCESS/PROCESS/process/models/tfcoil/base.py:255) dr_tf_plasma_case too small to accommodate the WP, forced to minimum value
 
******************************************* End of PROCESS Output ********************************************
 

Plot single solution¶

First we will look at the original large tokamak optimum solution. We will plot its solution, showing optimisation parameters normalised to their initial values.

In [4]:

large_tokamak_mfile = working_dir / "large_tokamak_MFILE.DAT"
# Plot the solution
runs_metadata = [
    RunMetadata(large_tokamak_mfile, "large tokamak"),
]

# Figure and dataframe returned for optional further modification
fig1, df1 = plot_mfile_solutions(
    runs_metadata=runs_metadata,
    plot_title="Large tokamak solution",
)
df1

large_tokamak_mfile = working_dir / "large_tokamak_MFILE.DAT" # Plot the solution runs_metadata = [ RunMetadata(large_tokamak_mfile, "large tokamak"), ] # Figure and dataframe returned for optional further modification fig1, df1 = plot_mfile_solutions( runs_metadata=runs_metadata, plot_title="Large tokamak solution", ) df1

Out[4]:

	tag	minmax	objf_name	norm_objf	itvar001_name	xcm001	itvar002_name	xcm002	itvar003_name	xcm003	...	itvar015_name	xcm015	itvar016_name	xcm016	itvar017_name	xcm017	itvar018_name	xcm018	itvar019_name	xcm019
0	large tokamak	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.878188	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.321601	...	c_tf_turn	1.384599	f_nd_alpha_electron	0.670216	f_a_cs_turn_steel	0.902271	f_nd_impurity_electrons(13)	3.022193	dr_tf_wp_with_insulation	1.051228

1 rows × 42 columns

Comparing optimum solutions¶

Now we will see the effect that varying an input parameter in the large tokamak input file has on the optimum solution found.

Here, the minimum allowable value for net electric power, p_plant_electric_net_required_mw, has been changed from 400MW to 200MW, and PROCESS has found a different optimum solution.

We can plot the two MFILEs together, showing normalised values of the optimisation parameters at the solution points, as well as the objective function values.

In [5]:

large_tokamak_varied_min_net_electric_mfile = (
    working_dir / "large_tokamak_varied_min_net_electric_MFILE.DAT"
)
runs_metadata = [
    RunMetadata(large_tokamak_mfile, "original"),
    RunMetadata(
        large_tokamak_varied_min_net_electric_mfile,
        "changed min net electric",
    ),
]

fig2, df2 = plot_mfile_solutions(
    runs_metadata=runs_metadata,
    plot_title="2 large tokamak solutions",
)
df2

large_tokamak_varied_min_net_electric_mfile = ( working_dir / "large_tokamak_varied_min_net_electric_MFILE.DAT" ) runs_metadata = [ RunMetadata(large_tokamak_mfile, "original"), RunMetadata( large_tokamak_varied_min_net_electric_mfile, "changed min net electric", ), ] fig2, df2 = plot_mfile_solutions( runs_metadata=runs_metadata, plot_title="2 large tokamak solutions", ) df2

Out[5]:

	tag	minmax	objf_name	norm_objf	itvar001_name	xcm001	itvar002_name	xcm002	itvar003_name	xcm003	...	itvar015_name	xcm015	itvar016_name	xcm016	itvar017_name	xcm017	itvar018_name	xcm018	itvar019_name	xcm019
0	original	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.878188	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.321601	...	c_tf_turn	1.384599	f_nd_alpha_electron	0.670216	f_a_cs_turn_steel	0.902271	f_nd_impurity_electrons(13)	3.022193	dr_tf_wp_with_insulation	1.051228
1	changed min net electric	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.837349	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.054697	...	c_tf_turn	1.190939	f_nd_alpha_electron	0.777458	f_a_cs_turn_steel	0.879007	f_nd_impurity_electrons(13)	0.978412	dr_tf_wp_with_insulation	1.128928

2 rows × 42 columns

We can compare the inequality constraint equations for both solutions.

In [6]:

import matplotlib.pyplot as plt

import process.core.io.mfile as mf
from process.core.io.plot_proc import plot_inequality_constraint_equations

original_mfile = mf.MFile((large_tokamak_mfile).as_posix())
new_mfile = mf.MFile((large_tokamak_varied_min_net_electric_mfile).as_posix())

f, axs = plt.subplots(1, 2)
axs[0].set_position([0.0, 0.0, 1.2, 1.5])
axs[1].set_position([1.9, 0.0, 1.2, 1.5])
plot_inequality_constraint_equations(axis=axs[0], m_file=original_mfile, scan=-1)
plot_inequality_constraint_equations(axis=axs[1], m_file=new_mfile, scan=-1)
axs[0].set_title("Original optimum solution")
axs[1].set_title("New optimum solution when changing min net electric")
f.suptitle("Inequality Constraint Equations", y=1.6, x=1.4)

import matplotlib.pyplot as plt import process.core.io.mfile as mf from process.core.io.plot_proc import plot_inequality_constraint_equations original_mfile = mf.MFile((large_tokamak_mfile).as_posix()) new_mfile = mf.MFile((large_tokamak_varied_min_net_electric_mfile).as_posix()) f, axs = plt.subplots(1, 2) axs[0].set_position([0.0, 0.0, 1.2, 1.5]) axs[1].set_position([1.9, 0.0, 1.2, 1.5]) plot_inequality_constraint_equations(axis=axs[0], m_file=original_mfile, scan=-1) plot_inequality_constraint_equations(axis=axs[1], m_file=new_mfile, scan=-1) axs[0].set_title("Original optimum solution") axs[1].set_title("New optimum solution when changing min net electric") f.suptitle("Inequality Constraint Equations", y=1.6, x=1.4)

Out[6]:

Text(1.4, 1.6, 'Inequality Constraint Equations')

To have lower net electric, PROCESS has found a solution where:

• the fusion power has dropped, therefore the neutron wall load has gone down
• the toroidal field required has slightly dropped, therefore the case stress limits are lower as less current in the coils is needed

Other solution comparison plots¶

There are some other ways that you can compare solutions in PROCESS. We will demonstrate these for the same two MFILEs as above, but you can add in more MFILEs if you want.

Here we refer to the original large tokamak file as Large tokamak 1, and the new solution obtained by varying p_plant_electric_net_required_mw as Large tokamak 2.

Plot one solution normalised to another¶

Normalised differences, relative to the a given solution, can also be plotted.

In [7]:

runs_metadata = [
    RunMetadata(large_tokamak_mfile, "large tokamak 1"),
    RunMetadata(
        large_tokamak_varied_min_net_electric_mfile,
        "large tokamak 2",
    ),
]

fig3, df3 = plot_mfile_solutions(
    runs_metadata=runs_metadata,
    plot_title="Large tokamak 2 solution, relative to large tokamak 1",
    normalising_tag="large tokamak 1",
)
df3

runs_metadata = [ RunMetadata(large_tokamak_mfile, "large tokamak 1"), RunMetadata( large_tokamak_varied_min_net_electric_mfile, "large tokamak 2", ), ] fig3, df3 = plot_mfile_solutions( runs_metadata=runs_metadata, plot_title="Large tokamak 2 solution, relative to large tokamak 1", normalising_tag="large tokamak 1", ) df3

/home/runner/work/PROCESS/PROCESS/process/core/io/plot_solutions.py:116: UserWarning: Double-normalising: using opt params normalised to each solution and normalising again to another solution. Are you sure?
  warn(

Out[7]:

	tag	minmax	objf_name	norm_objf	itvar001_name	xcm001	itvar002_name	xcm002	itvar003_name	xcm003	...	itvar015_name	xcm015	itvar016_name	xcm016	itvar017_name	xcm017	itvar018_name	xcm018	itvar019_name	xcm019
0	large tokamak 1	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.878188	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.321601	...	c_tf_turn	1.384599	f_nd_alpha_electron	0.670216	f_a_cs_turn_steel	0.902271	f_nd_impurity_electrons(13)	3.022193	dr_tf_wp_with_insulation	1.051228
1	large tokamak 2	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.837349	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.054697	...	c_tf_turn	1.190939	f_nd_alpha_electron	0.777458	f_a_cs_turn_steel	0.879007	f_nd_impurity_electrons(13)	0.978412	dr_tf_wp_with_insulation	1.128928

2 rows × 42 columns

RMS Errors¶

Plot RMS errors of multiple solutions relative to a reference solution.

In [8]:

fig5, df5 = plot_mfile_solutions(
    runs_metadata,
    "Large tokamak 2 solution with RMS errors normalised to large tokamak 1",
    normalising_tag="large tokamak 1",
    rmse=True,
)
df5

fig5, df5 = plot_mfile_solutions( runs_metadata, "Large tokamak 2 solution with RMS errors normalised to large tokamak 1", normalising_tag="large tokamak 1", rmse=True, ) df5

/home/runner/work/PROCESS/PROCESS/process/core/io/plot_solutions.py:116: UserWarning: Double-normalising: using opt params normalised to each solution and normalising again to another solution. Are you sure?
  warn(

Out[8]:

	tag	minmax	objf_name	norm_objf	itvar001_name	xcm001	itvar002_name	xcm002	itvar003_name	xcm003	...	itvar015_name	xcm015	itvar016_name	xcm016	itvar017_name	xcm017	itvar018_name	xcm018	itvar019_name	xcm019
0	large tokamak 1	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.878188	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.321601	...	c_tf_turn	1.384599	f_nd_alpha_electron	0.670216	f_a_cs_turn_steel	0.902271	f_nd_impurity_electrons(13)	3.022193	dr_tf_wp_with_insulation	1.051228
1	large tokamak 2	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.837349	rmajor	1.0	temp_plasma_electron_vol_avg_kev	1.054697	...	c_tf_turn	1.190939	f_nd_alpha_electron	0.777458	f_a_cs_turn_steel	0.879007	f_nd_impurity_electrons(13)	0.978412	dr_tf_wp_with_insulation	1.128928

2 rows × 42 columns

Solutions normalised by range¶

Use nitvar values instead; the solution optimisation parameters are normalised to the range of their upper and lower bounds.

In [9]:

fig6, df6 = plot_mfile_solutions(
    runs_metadata,
    "Large tokamak 2 solution normalised to the range of the optimisation parameters",
    normalisation_type="range",
)
df6

fig6, df6 = plot_mfile_solutions( runs_metadata, "Large tokamak 2 solution normalised to the range of the optimisation parameters", normalisation_type="range", ) df6

Out[9]:

	tag	minmax	objf_name	norm_objf	itvar001_name	nitvar001	itvar002_name	nitvar002	itvar003_name	nitvar003	...	itvar015_name	nitvar015	itvar016_name	nitvar016	itvar017_name	nitvar017	itvar018_name	nitvar018	itvar019_name	nitvar019
0	large tokamak 1	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.166578	rmajor	3.376943e-10	temp_plasma_electron_vol_avg_kev	0.109569	...	c_tf_turn	0.999956	f_nd_alpha_electron	0.340432	f_a_cs_turn_steel	0.759554	f_nd_impurity_electrons(13)	0.114842	dr_tf_wp_with_insulation	0.078509
1	large tokamak 2	1.0	major radius	1.6	b_plasma_toroidal_on_axis	0.158816	rmajor	2.663665e-10	temp_plasma_electron_vol_avg_kev	0.075675	...	c_tf_turn	0.496441	f_nd_alpha_electron	0.554915	f_a_cs_turn_steel	0.739943	f_nd_impurity_electrons(13)	0.037179	dr_tf_wp_with_insulation	0.102790

2 rows × 42 columns

Actual values¶

In [10]:

fig7, df7 = plot_mfile_solutions(
    runs_metadata,
    "Actual values of optimisation parameters for large tokamak 1 and 2 solutions",
    normalisation_type=None,
)
df7

fig7, df7 = plot_mfile_solutions( runs_metadata, "Actual values of optimisation parameters for large tokamak 1 and 2 solutions", normalisation_type=None, ) df7

Out[10]:

	tag	minmax	objf_name	norm_objf	itvar001_name	itvar001	itvar002_name	itvar002	itvar003_name	itvar003	...	itvar015_name	itvar015	itvar016_name	itvar016	itvar017_name	itvar017	itvar018_name	itvar018	itvar019_name	itvar019
0	large tokamak 1	1.0	major radius	1.6	b_plasma_toroidal_on_axis	5.005673	rmajor	8.0	temp_plasma_electron_vol_avg_kev	15.859208	...	c_tf_turn	89998.902736	f_nd_alpha_electron	0.067022	f_a_cs_turn_steel	0.721817	f_nd_impurity_electrons(13)	0.001148	dr_tf_wp_with_insulation	0.525614
1	large tokamak 2	1.0	major radius	1.6	b_plasma_toroidal_on_axis	4.772888	rmajor	8.0	temp_plasma_electron_vol_avg_kev	12.656366	...	c_tf_turn	77411.013813	f_nd_alpha_electron	0.077746	f_a_cs_turn_steel	0.703206	f_nd_impurity_electrons(13)	0.000372	dr_tf_wp_with_insulation	0.564464

2 rows × 42 columns

In [11]:

temp_dir.cleanup()

temp_dir.cleanup()