Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions documentation/source/development/add-vars.md
Original file line number Diff line number Diff line change
Expand Up @@ -142,13 +142,13 @@ After following the instruction to add an input variable, you can make the varia

## Add a constraint equation

Constraint equations are added to *PROCESS* in the `process/constraints.py` file. They are registered with the `ConstraintManager` whenever the application is run. Each equation has a unique name that is currently an integer, however upgrades to the input file format in the future will allow arbitrary hashable constraint names.
Constraint equations are added to *PROCESS* in the `process/core/optimisation/constraints.py` file. They are registered with the `ConstraintManager` whenever the application is run. Each equation has a unique name that is currently an integer, however upgrades to the input file format in the future will allow arbitrary hashable constraint names.

A constraint is simply added by registering the constraint to the manager using a decorator.

```python
@ConstraintManager.register_constraint(1234, "m", "=")
def my_constraint_function(): ...
def my_constraint_function(constraint_registration): ...
```
The arguments to the `register_constraint` function are:

Expand All @@ -161,13 +161,13 @@ The arguments to the `register_constraint` function are:

- Normalised residual error
- Constraint value
- Constraint error
- Constraint bound
- Constraint residual

The recommended way to do this is using one of the functions `geq`, `leq`, or `eq` depending on whether the constraint is desired to be $v\geq b$, $v\leq b$, or $v=b$, respectively.

```python
@ConstraintManager.register_constraint(1234, "m", "=")
def my_constraint_function():
normalised_residual = ...
value = ...
error = ...
return ConstraintResult(normalised_residual, value, error)
def my_constraint_function(constraint_registration):
return geq(value, bound, constraint_registration)
```
2 changes: 1 addition & 1 deletion documentation/source/io/input-guide.md
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ If no input file name is given as the first argument to the code, it assumes an
## Constraints

PROCESS permits a large number of constraint equations, all of which are formulated
in the source file `constraint equations.f90`.
in the source file `process/core/optimisation/constraints.py`.

In the input file constraint equations are specified as in the following example:

Expand Down
75 changes: 37 additions & 38 deletions documentation/source/solver/solver-guide.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,16 +3,40 @@

Any computer program naturally contains many equations. The built-in equation
solvers within PROCESS act on a special class, known as constraint equations,
all of which are formulated in the source file `constraint equations.f90`. These
all of which are formulated in the source file `process/core/optimisation/constraints.py`. These
can be split into two types:

**Equality constraints (consistency equations)** -- that enforce consistency between the physics and
engineering parameters in the various models
**Equality constraints (consistency equations)** that enforce consistency between the physics and
engineering parameters in the various models.

**Inequality constraints (limit equations)** -- that enforce various parameters to lie within their allowed
limits. The `neqns` constraint equations that the user chooses for a given run are
activated by including the equation numbers in the first `neqns` elements of
array `icc`.
**Inequality constraints (limit equations)** that enforce a value be greater/less than or equal to some bound (limit).

A PROCESS input file will, for example, define which constraint equations are being used as follows:

```
...
neqns = 3

* Equalities
icc = 1 * Beta
icc = 2 * Global power balance
icc = 11 * Radial build

* Inequalities
icc = 9 * Fusion power upper limit
icc = 5 * Density upper limit
icc = 24 * Beta upper limit
icc = 15 * LH power threshold limit
...
```

Here each `icc=n` statement tells PROCESS to activate a constraint with the name `n`. A list of the constraints and
their corresponding names can be found [here](../../source/reference/process/data_structure/numerics/#process.data_structure.numerics.lablcc).

The `neqns = 3` statement is telling PROCESS to treat the first `3` equations as equality constraints, and the rest as inequality constraints. Therefore, it is imperative that all equality constraints are stated before any inequality constraints.


In both types of equations, an optimiser/solver uses the normalised residuals $c_i$ of the constraints (and sometimes its gradient, depending on the solver/optimiser) to guide the solution towards one that satisfies all of the constraints.

## Consistency Equations

Expand All @@ -30,59 +54,34 @@ $$
c_i = 1 - \frac{g}{h}
$$

The equation solvers VMCON and HYBRD need the constraint equations $c_i$ to be given in
the form shown, since they adjust the iteration variables so as to obtain $c_i = 0$,
thereby ensuring that $g = h$.
The optimiser/solver will attempt to find a solution that produces $c_i = 0$ for all equality constraints.

## Limit Equations

The limit equations are inequalities that ensure that various physics or engineering
limits are not exceeded.

The f-values are used as follows. In general, limit equations have the form

$$
\mathrm{calculated\ quantity} = f \times \mathrm{maximum\ allowable\ value}
$$

where $f$ is the `f-value`. In optimisation mode, all iteration variables have prescribed
lower and upper bounds. If $f$ is chosen to be an iteration variable and is given a
lower bound of zero and an upper bound of one, then the limit equation does indeed
constrain the calculated quantity to lie between zero and its maximum allowable value,
as required.

As with the consistency equations, the general form of the limit equations is

$$
c_i = 1 - f.\frac{h_{max}}{h}
c_i = 1 - \frac{h_{max}}{h}
$$

where $h_{max}$ is the maximum allowed value of the quantity $h$. Sometimes, the limit
equation and `f-value` are used to ensure that quantity $h$ is larger than its minimum
value $h_{min}$. In this case, $0 ≤ f ≤ 1$ (as before), but the equation takes the form
where $h_{max}$ is the maximum allowed value of the quantity $h$, or

$$
c_i = 1 - f.\frac{h}{h_{min}}
$$

By fixing the `f-value` (i.e. not including it in the `ixc` array), the limit equations
can be used as equality constraints.
where $h_{min}$ is the minimum allowed value of the quantity $h$.

For example, to set the net electric power to a certain value, the following
should be carried out:

1. Activate `constraint 16` (net electric power lower limit) by including it in the `icc` array
2. Set `p_plant_electric_net_required_mw` to the required net electric power.

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.
The optimiser/solver will attempt to find a solution that produces $c_i \geq 0$ for all inequality constraints.

## Iteration Variables

Expand Down Expand Up @@ -111,7 +110,7 @@ Switch `ioptimz` should be set to 1 for optimisation mode.

If `ioptimz = 0`, a non-optimisation pass is performed first. Occasionally this provides a feasible set of initial conditions that aids convergence of the optimiser, but it is recommended to use `ioptimz = 1`.

Enable all the relevant consistency equations, and it is advisable to enable the corresponding iterations variables. A number of limit equations (inequality constraints) can also be activated. For limit equations, the corresponding f-value must be selected as an iteration variable. In optimisation more, the number of iteration variables is unlimited.
Enable all the relevant consistency equations, and it is advisable to enable the corresponding iterations variables. A number of limit equations (inequality constraints) can also be activated. In optimisation mode, the number of iteration variables is unlimited.

It may still be difficult, if not impossible, to reconcile the fusion power and the net electric power with the required values. This may well be due to the power conversion efficiency values being used.

Expand Down
2 changes: 1 addition & 1 deletion process/core/scan.py
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ def _missing_(cls, var):
"f_c_plasma_bootstrap_max", "Bootstrap_fraction", 12
)
boundu10 = ScanVariable("boundu(10)", "H_factor_upper_bound", 13)
fiooic = ScanVariable("fiooic", "TFC_Iop_/_Icrit_f-value", 14)
fiooic = ScanVariable("fiooic", "TFC_Iop_/_Icrit_margin", 14)
rmajor = ScanVariable("rmajor", "Plasma_major_radius_(m)", 16)
b_tf_inboard_max = ScanVariable("b_tf_inboard_max", "Max_toroidal_field_(T)", 17)
eta_cd_norm_hcd_primary_max = ScanVariable(
Expand Down
6 changes: 3 additions & 3 deletions tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT
Original file line number Diff line number Diff line change
Expand Up @@ -186,7 +186,7 @@ nd_plasma_electrons_vol_avg = 6.88360041658364314e+19
*ixc = 9
fdene = 1.2
*boundu(9) = 1.2
* DESCRIPTION: f-value for density limit (used to set max greenwald fraction)
* DESCRIPTION: margin for density limit (used to set max greenwald fraction)
* JUSTIFICATION: Used with icc=5 to enforce density limit

*ixc = 11 * pheat
Expand Down Expand Up @@ -242,7 +242,7 @@ f_c_plasma_non_inductive = 4.92677875380845121e-01
* JUSTIFICATION: We have a pulsed reactor so this can vary

fiooic = 1.0
* DESCRIPTION: f-value for TF coil operating current / critical current density ratio
* DESCRIPTION: margin for TF coil operating current / critical current density ratio
* JUSTIFICATION: Constraint equation 33 is used

ixc = 56 * tdmptf
Expand Down Expand Up @@ -282,7 +282,7 @@ dr_shld_vv_gap_inboard = 0.02
f_h_mode_margin = 1.1
*boundl(103) = 1.1
*boundu(103) = 1.2
* DESCRIPTION: f-value for L-H Power Threshold (now used to set lower bound)
* DESCRIPTION: margin for L-H Power Threshold (now used to set lower bound)
* JUSTIFICATION: Used with icc=15 to enforce H-mode

ixc = 109
Expand Down
Loading