PROCESS scan functionality behaves unexpectedly

[CoronelBuendia]

I have been performing a series of runs in PROCESS, for a typical conventional aspect ratio large tokamak problem, minimising R_0 for A in 2.6..3.6.

I used to do these kinds of runs frequently, prior to PROCESS being open-sourced. I would use the scan option liberally, usually just in one dimension. I found this feature useful, primarily because the run-time was better (the previous solution in a scan is usually a good starting point for a subsequent point).

Often, the scan would fail at a given point, and I would just start a new scan with a different starting point and "approach" that failure from a different direction.

So I took it as a given when now starting to do such scans in PROCESS v3.2.1 that I should be scanning. However this fails (using --solver vmcon_bounded), even for moderate deviations from the starting solution.

In fact, for this particular scan, I can find solutions for 100 points using the same starting vector (which I'm not sure would have worked in the last version of PROCESS I was using), but scanning I can't get much more than 12 or 13 points consecutively. I tried a range of different starting points and got these errors:

(/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/tf_coil.py:246) dr_tf_plasma_case to small to accommodate the WP, forced to minimum value
Starting scan point 13 of 100 : Aspect_ratio , aspect = 3.03939394
/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/problems/problem.py:1539: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  warnings.warn(
Traceback (most recent call last):-04
  File "/home/matti/miniforge/envs/bluemira/bin/process", line 7, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 753, in main
    Process(args)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 127, in __init__
    self.run_mode()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 231, in run_mode
    self.run.run()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 390, in run
    self.run_scan()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 474, in run_scan
    self.scan = Scan(self.models, self.solver)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 161, in __init__
    self.run_scan()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 190, in run_scan
    self.scan_1d()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 704, in scan_1d
    ifail = self.doopt()
            ^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 194, in doopt
    ifail = self.solver_handler.run()
            ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/solver_handler.py", line 50, in run
    ifail = self.solver.solve()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/solver.py", line 213, in solve
    x, _, _, res = solve(
                   ^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/pyvmcon/vmcon.py", line 172, in solve
    delta, lamda_equality, lamda_inequality = solve_qsp(
                                              ^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/pyvmcon/vmcon.py", line 331, in solve_qsp
    qsp.solve(**{"solver": cp.OSQP, **options})
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/problems/problem.py", line 609, in solve
    return solve_func(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/problems/problem.py", line 1191, in _solve
    data, solving_chain, inverse_data = self.get_problem_data(
                                        ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/problems/problem.py", line 755, in get_problem_data
    solving_chain = self._construct_chain(
                    ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/problems/problem.py", line 1009, in _construct_chain
    return construct_solving_chain(self, candidate_solvers, gp=gp,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/reductions/solvers/solving_chain.py", line 232, in construct_solving_chain
    reductions = _reductions_for_problem_class(problem, candidates, gp, solver_opts)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/cvxpy/reductions/solvers/solving_chain.py", line 147, in _reductions_for_problem_class
    raise DCPError(
cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:
QuadForm(var13357, [[7914.91 6591.39 ... 1.92 0.00]
 [6591.39 5498.00 ... 1.69 0.00]
 ...
 [1.92 1.69 ... 0.88 0.00]
 [0.00 0.00 ... 0.00 1.00]])

and

process.solver - CRITICAL - QSP failed to solve, indicating no feasible solution could be found.

 
********************************* PROCESS COULD NOT FIND A FEASIBLE SOLUTION *********************************
 
 
process.scan - CRITICAL - Solver returns with ifail /= 1. ifail=5
 
Traceback (most recent call last):
  File "/home/matti/miniforge/envs/bluemira/bin/process", line 7, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 753, in main
    Process(args)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 127, in __init__
    self.run_mode()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 231, in run_mode
    self.run.run()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 390, in run
    self.run_scan()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/main.py", line 474, in run_scan
    self.scan = Scan(self.models, self.solver)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 161, in __init__
    self.run_scan()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 190, in run_scan
    self.scan_1d()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/scan.py", line 706, in scan_1d
    write_output_files(models=self.models, ifail=ifail)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/caller.py", line 368, in write_output_files
    caller.call_models_and_write_output(
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/caller.py", line 137, in call_models_and_write_output
    finalise(self.models, ifail)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/final.py", line 40, in finalise
    op.write(models, constants.NOUT)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/output.py", line 44, in write
    models.physics.outplas()
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/physics.py", line 6102, in outplas
    self.output_confinement_comparison(istell=stellarator_variables.istell)
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/physics.py", line 6628, in output_confinement_comparison
    physics_variables.hfac[i_confinement_time - 1] = self.find_other_h_factors(
                                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/process/physics.py", line 7582, in find_other_h_factors
    return root_scalar(fhz, bracket=(0.01, 150), xtol=0.001).root
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/scipy/optimize/_root_scalar.py", line 286, in root_scalar
    r, sol = methodc(f, a, b, args=args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/matti/miniforge/envs/bluemira/lib/python3.11/site-packages/scipy/optimize/_zeros_py.py", line 798, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

@timothy-nunn has looked into this a bit already and found that the Hessian approximation in PyVMCON is apparently non positive semi-definite.

Generally speaking, this doesn't make any sense to me: starting a local optimiser closer to a local optimum is almost guaranteed to speed up convergence... Thought I would flag this, as it tripped me up - big time. I can pass whoever is assigned this the IN.DAT I was using.

[timothy-nunn]

A couple of ideas that I have to investigate that might shed some light on this problem:

1. Perform a sweep with sweep = x, x, x, ..., x. One would expect that each scan point would converge to the same/very similar optima. In practice, it could be that starting so close to the optima--a scan starts at the previous optimal x--causes a drift over time which eventually leads the solver into a locally unsolvable area (QSP solver fails).
2. Try and perform this sweep with the old VMCON to see if this is an issue with PyVMCON or an issue with PROCESS (e.g. bad model behaviour).
3. Check how the hessian approximation matrix in PyVMCON is becoming non-positive semi-definite. There is a provision in PyVMCON to stop this happening, so is this provision not guaranteed or is this caused by numeric instability (e.g. floating point rounding error).

[CoronelBuendia]

I tried (1) for 15 points and it behaves as expected there is very little/no drift in solution vector over the scan points. All points after the first converge in 8 iterations (with the same convergence parameter to 3 SD as of the 3 rd point).

[timothy-nunn]

@chris-ashe please mark your PR for #3674 to solve this issue. I think the improvements to PyVMCON (and how PROCESS uses it) + the removal of two common crashing modes constitutes enough of an improvement

[timothy-nunn]

I have marked #4018 as closing this issue. As stated above, once this PR is merged, we will have fixed multiple problems observed when scanning:

• Propagate and then removed the cause of the DCPError (fixed in PyVMCON)
• Change of PyVMCON inner solver and tolerance to improve the quality of steps taken by the optimiser
• Removal of crash mode when H-factor could not be calculated (this was for output only and so it was sufficient to catch the error and set the variable to NaN)
• Removal of another crash mode when a temperature margin cannot be calculated.