Advanced QS Optimization
In this tutorial we will show an example of “precise” QS optimization using a multigrid approach, and using constrained optimization
[1]:
import sys
import os
sys.path.insert(0, os.path.abspath("."))
sys.path.append(os.path.abspath("../../../"))
If you have access to a GPU, uncomment the following two lines before any DESC or JAX related imports. You should see about an order of magnitude speed improvement with only these two lines of code!
[2]:
# from desc import set_device
# set_device("gpu")
As mentioned in DESC Documentation on performance tips, one can use compilation cache directory to reduce the compilation overhead time. Note: One needs to create jax-caches folder manually.
[ ]:
# import jax
# jax.config.update("jax_compilation_cache_dir", "../jax-caches")
# jax.config.update("jax_persistent_cache_min_entry_size_bytes", -1)
# jax.config.update("jax_persistent_cache_min_compile_time_secs", 0)
[3]:
import numpy as np
from desc.continuation import solve_continuation_automatic
from desc.equilibrium import EquilibriaFamily, Equilibrium
from desc.geometry import FourierRZToroidalSurface
from desc.grid import LinearGrid
from desc.objectives import (
AspectRatio,
FixBoundaryR,
FixBoundaryZ,
FixCurrent,
FixPressure,
FixPsi,
ForceBalance,
ObjectiveFunction,
QuasisymmetryTwoTerm,
)
from desc.optimize import Optimizer
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
DESC version 0.13.0+702.ge6b8a02dc.dirty,using JAX backend, jax version=0.4.33, jaxlib version=0.4.33, dtype=float64
Using device: CPU, with 56.70 GB available memory
Initial Guess
We start by creating an initial equilibrium and solving a standard fixed boundary problem:
[4]:
# create initial surface. Aspect ratio ~ 8, circular cross section with slight
# axis torsion to make it nonplanar
surf = FourierRZToroidalSurface(
R_lmn=[1, 0.125, 0.1],
Z_lmn=[-0.125, -0.1],
modes_R=[[0, 0], [1, 0], [0, 1]],
modes_Z=[[-1, 0], [0, -1]],
NFP=4,
)
# create initial equilibrium. Psi chosen to give B ~ 1 T. Could also give profiles here,
# default is zero pressure and zero current
eq = Equilibrium(M=4, N=4, Psi=0.04, surface=surf)
# this is usually all you need to solve a fixed boundary equilibrium
eq0 = solve_continuation_automatic(eq, verbose=0)[-1]
# it will be helpful to store intermediate results
eqfam = EquilibriaFamily(eq0)
Multigrid method with proximal optimizer
By “multigrid” method we mean we will start by optimizing over boundary modes with \(|m|, |n| \leq 1\), then \(|m|, |n| \leq 2\) and so on. To do this we’ll define a helper function that will create the necessary constraints and objectives for a given maximum mode number \(k\).
By a “proximal” optimizer we mean one that handles the equilibrium constraint by re-solving a fixed boundary equilibrium problem at each step, given the current position of the boundary. This is made more efficient by using a perturbed estimate based on the previous step as a warm start to the equilibrium sub-problem.
[5]:
def run_qh_step(k, eq):
"""Run a step of the precise QH optimization example from Landreman & Paul."""
# this step will only optimize boundary modes with |m|,|n| <= k
# create grid where we want to minimize QS error. Here we do it on 3 surfaces
grid = LinearGrid(
M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, rho=np.array([0.6, 0.8, 1.0]), sym=True
)
# we create an ObjectiveFunction, in this case made up of multiple objectives
# which will be combined in a least squares sense
objective = ObjectiveFunction(
(
# pass in the grid we defined, and don't forget the target helicity!
QuasisymmetryTwoTerm(eq=eq, helicity=(1, eq.NFP), grid=grid),
# try to keep the aspect ratio about the same
AspectRatio(eq=eq, target=8, weight=100),
),
)
# as opposed to SIMSOPT and STELLOPT where variables are assumed fixed, in DESC
# we assume variables are free. Here we decide which ones to fix, starting with
# the major radius (R mode = [0,0,0]) and all modes with m,n > k
R_modes = np.vstack(
(
[0, 0, 0],
eq.surface.R_basis.modes[
np.max(np.abs(eq.surface.R_basis.modes), 1) > k, :
],
)
)
Z_modes = eq.surface.Z_basis.modes[
np.max(np.abs(eq.surface.Z_basis.modes), 1) > k, :
]
# next we create the constraints, using the mode number arrays just created
# if we didn't pass those in, it would fix all the modes (like for the profiles)
constraints = (
ForceBalance(eq=eq),
FixBoundaryR(eq=eq, modes=R_modes),
FixBoundaryZ(eq=eq, modes=Z_modes),
FixPressure(eq=eq),
FixCurrent(eq=eq),
FixPsi(eq=eq),
)
# this is the default optimizer, which re-solves the equilibrium at each step
optimizer = Optimizer("proximal-lsq-exact")
eq_new, history = eq.optimize(
objective=objective,
constraints=constraints,
optimizer=optimizer,
maxiter=20, # we don't need to solve to optimality at each multigrid step
verbose=3,
copy=True, # don't modify original, return a new optimized copy
options={
# Sometimes the default initial trust radius is too big, allowing the
# optimizer to take too large a step in a bad direction. If this happens,
# we can manually specify a smaller starting radius. Each optimizer has a
# number of different options that can be used to tune the performance.
# See the documentation for more info.
"initial_trust_ratio": 1.0,
},
)
return eq_new
Lets look at the initial field:
[6]:
from desc.plotting import plot_boozer_surface
plot_boozer_surface(eq0);
We see that it is vaguely QH like, which is why we’re targeting QH symmetry. Now let’s run the optimization in steps and look at the intermediate result after each step:
[7]:
eq1 = run_qh_step(1, eq0)
eqfam.append(eq1)
plot_boozer_surface(eq1);
Building objective: QS two-term
Precomputing transforms
Timer: Precomputing transforms = 557 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 22.9 ms
Timer: Objective build = 635 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 50.9 ms
Timer: Objective build = 60.9 ms
Timer: Proximal projection build = 1.27 sec
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 298 ms
Timer: Linear constraint projection build = 2.04 sec
Number of parameters: 8
Number of objectives: 460
Timer: Initializing the optimization = 3.67 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.946e+01 2.907e+00
1 5 2.765e+01 1.810e+00 1.726e-02 2.009e+00
2 6 2.493e+01 2.718e+00 6.421e-03 8.119e-01
3 7 2.222e+01 2.710e+00 1.592e-02 5.414e-01
4 8 1.925e+01 2.970e+00 3.283e-02 3.295e-01
5 9 1.789e+01 1.359e+00 4.940e-02 4.881e-01
6 10 1.784e+01 5.185e-02 4.410e-02 4.985e-01
7 11 1.713e+01 7.141e-01 1.576e-02 1.530e-01
8 12 1.669e+01 4.341e-01 9.818e-03 1.024e-01
9 13 1.634e+01 3.502e-01 6.521e-03 1.583e-01
10 14 1.582e+01 5.208e-01 1.367e-02 5.408e-01
11 15 1.504e+01 7.801e-01 2.132e-02 5.002e-01
12 17 1.431e+01 7.343e-01 1.097e-02 2.051e-01
13 18 1.353e+01 7.752e-01 1.267e-02 7.186e-01
14 19 1.341e+01 1.213e-01 2.185e-02 1.978e+00
15 20 1.132e+01 2.087e+00 1.879e-02 1.042e-01
16 21 1.055e+01 7.726e-01 8.481e-03 2.524e-01
17 22 9.683e+00 8.677e-01 1.600e-02 9.819e-01
18 23 8.376e+00 1.307e+00 1.614e-02 4.978e-01
19 25 7.859e+00 5.164e-01 9.695e-03 1.039e-01
20 26 7.324e+00 5.348e-01 1.653e-02 2.702e-01
Warning: Maximum number of iterations has been exceeded.
Current function value: 7.324e+00
Total delta_x: 1.441e-01
Iterations: 20
Function evaluations: 26
Jacobian evaluations: 21
Timer: Solution time = 50.9 sec
Timer: Avg time per step = 2.42 sec
==============================================================================================================
Start --> End
Total (sum of squares): 2.946e+01 --> 7.324e+00,
Maximum absolute Quasi-symmetry (1,4) two-term error: 6.229e-01 --> 4.966e-01 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.950e-04 --> 6.508e-04 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.501e-01 --> 7.247e-02 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 5.894e-01 --> 4.699e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.845e-04 --> 6.158e-04 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.420e-01 --> 6.857e-02 (normalized)
Aspect ratio: 8.000e+00 --> 7.998e+00 (dimensionless)
Maximum absolute Force error: 4.341e+01 --> 1.105e+03 (N)
Minimum absolute Force error: 4.127e-03 --> 7.947e-02 (N)
Average absolute Force error: 3.601e+00 --> 6.058e+01 (N)
Maximum absolute Force error: 4.261e-05 --> 1.085e-03 (normalized)
Minimum absolute Force error: 4.051e-09 --> 7.802e-08 (normalized)
Average absolute Force error: 3.535e-06 --> 5.948e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
==============================================================================================================
[8]:
eq2 = run_qh_step(2, eq1)
eqfam.append(eq2)
plot_boozer_surface(eq2);
Building objective: QS two-term
Precomputing transforms
Timer: Precomputing transforms = 48.4 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 22.1 ms
Timer: Objective build = 84.1 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 51.0 ms
Timer: Objective build = 61.2 ms
Timer: Proximal projection build = 437 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 289 ms
Timer: Linear constraint projection build = 1.43 sec
Number of parameters: 24
Number of objectives: 460
Timer: Initializing the optimization = 2.21 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.670e+01 1.087e+00
1 5 2.768e+01 9.015e+00 1.425e-02 1.742e+00
2 6 2.230e+01 5.378e+00 1.463e-02 2.167e+00
3 7 1.499e+01 7.318e+00 1.365e-02 1.287e+00
4 8 7.989e+00 6.998e+00 1.962e-02 1.378e+00
5 9 3.545e+00 4.444e+00 2.534e-02 1.145e+00
6 11 1.917e+00 1.627e+00 1.476e-02 4.953e-01
7 12 1.369e+00 5.487e-01 1.660e-02 7.344e-01
8 13 1.055e+00 3.140e-01 1.627e-02 5.301e-01
9 14 8.836e-01 1.710e-01 1.624e-02 4.852e-01
10 16 6.508e-01 2.328e-01 5.191e-03 7.269e-02
11 18 6.318e-01 1.909e-02 2.614e-03 2.247e-02
12 19 6.109e-01 2.083e-02 5.859e-03 1.140e-01
13 21 5.886e-01 2.228e-02 2.349e-03 2.397e-02
14 22 5.713e-01 1.735e-02 6.001e-03 1.089e-01
15 24 5.537e-01 1.763e-02 2.301e-03 2.033e-02
16 25 5.446e-01 9.016e-03 6.195e-03 9.574e-02
17 27 5.347e-01 9.994e-03 2.659e-03 1.620e-02
18 28 5.339e-01 7.206e-04 6.496e-03 3.371e-02
19 29 5.318e-01 2.142e-03 1.212e-03 1.422e-02
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 5.318e-01
Total delta_x: 1.144e-01
Iterations: 19
Function evaluations: 29
Jacobian evaluations: 20
Timer: Solution time = 35.3 sec
Timer: Avg time per step = 1.77 sec
==============================================================================================================
Start --> End
Total (sum of squares): 3.670e+01 --> 5.318e-01,
Maximum absolute Quasi-symmetry (1,4) two-term error: 4.966e-01 --> 9.640e-02 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 6.508e-04 --> 1.578e-05 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 7.247e-02 --> 1.005e-02 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.053e+00 --> 2.044e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.380e-03 --> 3.345e-05 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.536e-01 --> 2.131e-02 (normalized)
Aspect ratio: 7.998e+00 --> 8.000e+00 (dimensionless)
Maximum absolute Force error: 1.105e+03 --> 8.833e+02 (N)
Minimum absolute Force error: 7.947e-02 --> 1.136e-01 (N)
Average absolute Force error: 6.058e+01 --> 7.925e+01 (N)
Maximum absolute Force error: 1.624e-03 --> 1.298e-03 (normalized)
Minimum absolute Force error: 1.168e-07 --> 1.669e-07 (normalized)
Average absolute Force error: 8.903e-05 --> 1.165e-04 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
==============================================================================================================
[9]:
eq3 = run_qh_step(3, eq2)
eqfam.append(eq3)
plot_boozer_surface(eq3);
Building objective: QS two-term
Precomputing transforms
Timer: Precomputing transforms = 49.7 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 21.9 ms
Timer: Objective build = 82.9 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 50.7 ms
Timer: Objective build = 60.0 ms
Timer: Proximal projection build = 435 ms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 201 ms
Timer: Linear constraint projection build = 1.37 sec
Number of parameters: 48
Number of objectives: 460
Timer: Initializing the optimization = 2.07 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 4.410e-01 1.645e-01
1 4 3.370e-01 1.040e-01 8.384e-03 4.815e-01
2 5 2.438e-01 9.320e-02 4.402e-03 1.833e-01
3 6 2.236e-01 2.013e-02 5.653e-03 3.771e-01
4 7 1.406e-01 8.306e-02 1.699e-03 4.500e-02
5 8 1.397e-01 8.711e-04 3.333e-03 1.242e-01
6 9 1.201e-01 1.963e-02 9.391e-04 1.411e-02
7 10 1.149e-01 5.106e-03 1.757e-03 3.653e-02
8 11 1.149e-01 2.321e-05 1.674e-03 3.896e-02
9 12 1.118e-01 3.144e-03 4.968e-04 9.711e-03
10 13 1.103e-01 1.450e-03 4.784e-04 6.444e-03
11 14 1.096e-01 7.663e-04 4.539e-04 5.976e-03
12 40 1.096e-01 0.000e+00 0.000e+00 5.976e-03
Warning: A bad approximation caused failure to predict improvement.
Current function value: 1.096e-01
Total delta_x: 1.668e-02
Iterations: 12
Function evaluations: 40
Jacobian evaluations: 12
Timer: Solution time = 33.0 sec
Timer: Avg time per step = 2.54 sec
==============================================================================================================
Start --> End
Total (sum of squares): 4.410e-01 --> 1.096e-01,
Maximum absolute Quasi-symmetry (1,4) two-term error: 9.640e-02 --> 7.544e-02 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.578e-05 --> 3.806e-06 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.005e-02 --> 4.756e-03 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.861e-01 --> 1.456e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 3.046e-05 --> 7.347e-06 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.940e-02 --> 9.182e-03 (normalized)
Aspect ratio: 8.000e+00 --> 8.000e+00 (dimensionless)
Maximum absolute Force error: 8.833e+02 --> 6.867e+02 (N)
Minimum absolute Force error: 1.136e-01 --> 1.617e-02 (N)
Average absolute Force error: 7.925e+01 --> 5.656e+01 (N)
Maximum absolute Force error: 1.239e-03 --> 9.629e-04 (normalized)
Minimum absolute Force error: 1.593e-07 --> 2.267e-08 (normalized)
Average absolute Force error: 1.111e-04 --> 7.931e-05 (normalized)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
==============================================================================================================
We see that after only 3 multigrid steps we have achieved very straight contours of magnetic field strength. These could be further refined by running for more iterations, using higher resolution, tighter tolerances, etc.
As a final comparison, we’ll look at the maximum symmetry breaking boozer harmonic for each step of the equilibrium
[10]:
import matplotlib.pyplot as plt
from desc.plotting import plot_boozer_modes, plot_boundaries
fig, ax = plt.subplots()
colors = ["r", "g", "c", "m"]
for i, (eq, color) in enumerate(zip(eqfam, colors)):
plot_boozer_modes(
eq, color=color, helicity=(1, eq.NFP), max_only=True, label=f"Step {i}", ax=ax
);
Constrained Optimization
Next, we’ll do a similar optimization but this time treating it as a constrained optimization problem, where we attempt to minimize QS error subject to more complicated constraints. We’ll start with the same QS objective:
[11]:
# create grid where we want to minimize QS error. Here we do it on 3 surfaces
grid = LinearGrid(
M=eq0.M_grid, N=eq0.N_grid, NFP=eq0.NFP, rho=np.array([0.6, 0.8, 1.0]), sym=True
)
objective = ObjectiveFunction(
(
# pass in the grid we defined, and don't forget the target helicity!
QuasisymmetryTwoTerm(eq=eq0, helicity=(1, eq.NFP), grid=grid),
),
)
For constraints, we’ll include the standard force balance to start. In the previous example, fixing certain boundary modes served as a form of regularization to prevent the solution from going into a bad local minimum. In this case however, instead of fixing a range of boundary modes we will only fix the \(R_{00}\) mode, and include constraints on aspect ratio, volume, and elongation to keep the solution from going off in a bad direction.
Finally, we also include a constraint on the average rotational transform:
[12]:
from desc.objectives import Elongation, RotationalTransform, Volume
constraints = (
ForceBalance(eq=eq0),
# try to keep the aspect ratio between 7 and 9
AspectRatio(eq=eq0, bounds=(7, 9)),
# similarly, try to keep it from getting too elongated
Elongation(eq=eq0, bounds=(0, 3)),
# Keep volume the same as the initial volume
Volume(eq=eq0, target=eq0.compute("V")["V"]),
# target for average iota
RotationalTransform(eq=eq0, target=1.1, loss_function="mean"),
# fix major radius
FixBoundaryR(eq=eq0, modes=[0, 0, 0]),
# fix vacuum profiles
FixPressure(eq=eq0),
FixCurrent(eq=eq0),
FixPsi(eq=eq0),
)
Finally, we’ll use an optimizer that can handle general nonlinear constraints (the proximal-lsq-exact optimizer can only handle equilibrium constraints such as ForceBalance and regular linear constraints like Fix*). In this case we use a least-squares augmented Lagrangian method.
[13]:
optimizer = Optimizer("lsq-auglag")
eqa, history = eq0.optimize(
objective=objective,
constraints=constraints,
optimizer=optimizer,
# each iteration of the augmented Lagrangian optimizer is cheaper than a step of a
# proximal optimizer, but it generally requires more iterations to converge
maxiter=200,
copy=True,
verbose=3,
options={},
)
Building objective: QS two-term
Precomputing transforms
Timer: Precomputing transforms = 50.5 ms
Timer: Objective build = 58.9 ms
Building objective: lcfs R
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Timer: Objective build = 92.5 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 51.4 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 22.3 ms
Building objective: elongation
Precomputing transforms
Timer: Precomputing transforms = 22.4 ms
Building objective: volume
Precomputing transforms
Timer: Precomputing transforms = 22.6 ms
Building objective: rotational transform
Precomputing transforms
Timer: Precomputing transforms = 601 ms
Timer: Objective build = 846 ms
Timer: Linear constraint projection build = 2.90 sec
Timer: Linear constraint projection build = 95.0 ms
Number of parameters: 200
Number of objectives: 459
Number of equality constraints: 852
Number of inequality constraints: 2
Timer: Initializing the optimization = 4.21 sec
Starting optimization
Using method: lsq-auglag
Iteration Total nfev Cost Cost reduction Step norm Optimality Constr viol. Penalty param max(|mltplr|)
0 1 2.946e+01 7.110e+00 8.523e-01 1.000e+01 0.000e+00
1 2 2.674e+01 2.722e+00 4.197e-03 6.695e+00 8.531e-01 1.000e+01 0.000e+00
2 3 1.777e+01 8.965e+00 1.895e-02 5.197e+00 8.562e-01 1.000e+01 0.000e+00
3 4 1.878e+00 1.589e+01 8.647e-02 1.277e+00 8.566e-01 1.000e+01 0.000e+00
4 5 2.014e+00 -1.365e-01 3.885e-01 1.086e+00 5.454e-01 1.000e+01 0.000e+00
5 6 4.470e-01 1.567e+00 1.526e-01 1.176e-01 6.071e-01 1.000e+01 0.000e+00
6 7 3.192e-01 1.278e-01 1.291e-01 1.431e-01 5.978e-01 1.000e+01 0.000e+00
7 8 3.636e-01 -4.442e-02 3.051e-01 1.421e-01 4.838e-01 1.000e+01 0.000e+00
8 10 3.428e-01 2.083e-02 2.119e-01 2.889e-01 3.343e-01 1.000e+01 0.000e+00
9 12 4.590e-01 -1.162e-01 2.685e-01 4.732e-01 1.650e-01 1.000e+01 0.000e+00
10 14 6.385e-02 3.952e-01 8.674e-02 4.182e-02 1.402e-01 1.000e+01 0.000e+00
11 16 4.167e-02 2.218e-02 6.740e-02 2.541e-02 1.078e-01 1.000e+01 0.000e+00
12 18 3.656e-02 5.107e-03 7.832e-02 2.550e-02 8.356e-02 1.000e+01 0.000e+00
13 20 3.301e-02 3.553e-03 6.883e-02 2.508e-02 6.270e-02 1.000e+01 0.000e+00
14 22 3.014e-02 2.870e-03 6.000e-02 3.115e-02 4.825e-02 1.000e+01 0.000e+00
15 24 2.758e-02 2.557e-03 4.460e-02 2.829e-02 3.817e-02 1.000e+01 0.000e+00
16 26 2.530e-02 2.277e-03 3.299e-02 2.273e-02 3.146e-02 1.000e+01 0.000e+00
17 28 2.334e-02 1.964e-03 2.684e-02 2.272e-02 2.691e-02 1.000e+01 0.000e+00
18 30 2.157e-02 1.768e-03 2.437e-02 2.219e-02 2.357e-02 1.000e+01 0.000e+00
19 32 1.995e-02 1.622e-03 2.380e-02 2.121e-02 2.088e-02 1.000e+01 0.000e+00
20 34 1.848e-02 1.469e-03 2.433e-02 2.012e-02 1.859e-02 1.000e+01 0.000e+00
21 36 1.716e-02 1.323e-03 2.681e-02 1.929e-02 1.657e-02 1.000e+01 0.000e+00
22 38 1.597e-02 1.190e-03 2.715e-02 1.933e-02 1.474e-02 1.000e+01 0.000e+00
23 40 1.489e-02 1.080e-03 2.760e-02 2.038e-02 1.380e-02 1.000e+01 0.000e+00
24 42 1.391e-02 9.736e-04 2.799e-02 2.234e-02 1.329e-02 1.000e+01 0.000e+00
25 44 1.304e-02 8.716e-04 3.090e-02 2.460e-02 1.286e-02 1.000e+01 0.000e+00
26 46 1.227e-02 7.748e-04 3.065e-02 2.690e-02 1.270e-02 1.000e+01 0.000e+00
27 48 1.155e-02 7.215e-04 3.051e-02 2.831e-02 1.252e-02 1.000e+01 0.000e+00
28 50 1.085e-02 7.006e-04 3.000e-02 2.877e-02 1.232e-02 1.000e+01 0.000e+00
29 52 1.013e-02 7.151e-04 2.918e-02 2.817e-02 1.209e-02 1.000e+01 0.000e+00
30 54 9.379e-03 7.520e-04 2.815e-02 2.645e-02 1.183e-02 1.000e+01 0.000e+00
31 56 8.587e-03 7.913e-04 2.695e-02 2.369e-02 1.155e-02 1.000e+01 0.000e+00
32 58 7.777e-03 8.105e-04 2.563e-02 2.020e-02 1.150e-02 1.000e+01 0.000e+00
33 60 6.980e-03 7.967e-04 2.425e-02 1.645e-02 1.143e-02 1.000e+01 0.000e+00
34 62 6.228e-03 7.516e-04 2.288e-02 1.292e-02 1.132e-02 1.000e+01 0.000e+00
35 64 5.547e-03 6.811e-04 2.166e-02 1.002e-02 1.119e-02 1.000e+01 0.000e+00
36 66 4.958e-03 5.895e-04 2.067e-02 1.028e-02 1.102e-02 1.000e+01 0.000e+00
37 68 4.482e-03 4.757e-04 1.984e-02 1.353e-02 1.082e-02 1.000e+01 0.000e+00
38 70 4.133e-03 3.488e-04 1.899e-02 1.727e-02 1.057e-02 1.000e+01 0.000e+00
39 72 3.892e-03 2.415e-04 1.804e-02 2.030e-02 1.030e-02 1.000e+01 0.000e+00
40 74 3.710e-03 1.821e-04 1.699e-02 2.221e-02 1.001e-02 1.000e+01 0.000e+00
41 75 3.539e-03 1.704e-04 1.580e-02 2.247e-02 9.709e-03 1.000e+01 0.000e+00
42 76 3.344e-03 1.950e-04 1.449e-02 2.209e-02 9.412e-03 1.000e+01 0.000e+00
43 77 3.137e-03 2.072e-04 1.321e-02 2.215e-02 9.121e-03 1.000e+01 0.000e+00
44 78 2.937e-03 2.003e-04 1.210e-02 2.222e-02 8.837e-03 1.000e+01 0.000e+00
45 79 2.746e-03 1.912e-04 1.244e-02 2.200e-02 8.563e-03 1.000e+01 0.000e+00
46 80 2.569e-03 1.768e-04 1.194e-02 2.164e-02 8.301e-03 1.000e+01 0.000e+00
47 81 2.408e-03 1.607e-04 1.158e-02 2.122e-02 8.051e-03 1.000e+01 0.000e+00
48 82 2.263e-03 1.453e-04 1.134e-02 2.079e-02 7.812e-03 1.000e+01 0.000e+00
49 83 2.132e-03 1.312e-04 1.115e-02 2.034e-02 7.589e-03 1.000e+01 0.000e+00
50 84 2.013e-03 1.187e-04 1.097e-02 1.988e-02 7.394e-03 1.000e+01 0.000e+00
51 85 1.905e-03 1.083e-04 1.079e-02 1.937e-02 7.201e-03 1.000e+01 0.000e+00
52 86 1.804e-03 1.001e-04 1.061e-02 1.880e-02 7.012e-03 1.000e+01 0.000e+00
53 87 1.711e-03 9.335e-05 1.044e-02 1.819e-02 6.829e-03 1.000e+01 0.000e+00
54 88 1.624e-03 8.752e-05 1.026e-02 1.755e-02 6.652e-03 1.000e+01 0.000e+00
55 89 1.541e-03 8.271e-05 1.008e-02 1.690e-02 6.480e-03 1.000e+01 0.000e+00
56 90 1.461e-03 7.967e-05 9.895e-03 1.623e-02 6.311e-03 1.000e+01 0.000e+00
57 91 1.383e-03 7.864e-05 9.719e-03 1.549e-02 6.144e-03 1.000e+01 0.000e+00
58 92 1.304e-03 7.850e-05 9.562e-03 1.459e-02 5.977e-03 1.000e+01 0.000e+00
59 93 1.227e-03 7.714e-05 9.431e-03 1.349e-02 5.812e-03 1.000e+01 0.000e+00
60 95 1.154e-03 7.285e-05 9.326e-03 1.220e-02 5.649e-03 1.000e+01 0.000e+00
61 97 1.088e-03 6.564e-05 9.241e-03 1.083e-02 5.489e-03 1.000e+01 0.000e+00
62 99 1.031e-03 5.720e-05 9.168e-03 9.456e-03 5.331e-03 1.000e+01 0.000e+00
63 101 9.820e-04 4.928e-05 9.099e-03 8.159e-03 5.177e-03 1.000e+01 0.000e+00
64 103 9.389e-04 4.305e-05 9.027e-03 6.954e-03 5.029e-03 1.000e+01 0.000e+00
65 105 9.002e-04 3.877e-05 8.948e-03 6.285e-03 4.886e-03 1.000e+01 0.000e+00
66 107 8.641e-04 3.603e-05 8.862e-03 5.776e-03 4.751e-03 1.000e+01 0.000e+00
67 109 8.298e-04 3.431e-05 8.770e-03 5.255e-03 4.623e-03 1.000e+01 0.000e+00
68 111 7.965e-04 3.334e-05 8.672e-03 4.739e-03 4.505e-03 1.000e+01 0.000e+00
69 113 7.635e-04 3.300e-05 8.563e-03 4.242e-03 4.396e-03 1.000e+01 0.000e+00
70 115 7.307e-04 3.281e-05 8.452e-03 4.374e-03 4.298e-03 1.000e+01 0.000e+00
71 117 7.000e-04 3.068e-05 8.378e-03 5.088e-03 4.212e-03 1.000e+01 0.000e+00
72 119 6.784e-04 2.160e-05 8.423e-03 6.097e-03 4.138e-03 1.000e+01 0.000e+00
73 121 6.767e-04 1.678e-06 8.650e-03 7.271e-03 4.072e-03 1.000e+01 0.000e+00
74 123 6.809e-04 -4.197e-06 8.854e-03 8.027e-03 4.012e-03 1.000e+01 0.000e+00
75 125 6.475e-04 3.346e-05 8.677e-03 7.881e-03 3.959e-03 1.000e+01 0.000e+00
76 127 6.128e-04 3.468e-05 8.457e-03 7.273e-03 3.890e-03 1.000e+01 0.000e+00
77 129 6.026e-04 1.016e-05 8.596e-03 2.467e-02 3.808e-03 1.000e+01 0.000e+00
78 130 5.953e-04 7.287e-06 9.986e-03 1.006e-02 3.754e-03 1.000e+01 0.000e+00
79 131 5.828e-04 1.254e-05 1.043e-02 1.010e-02 3.707e-03 1.000e+01 0.000e+00
80 132 5.630e-04 1.983e-05 1.103e-02 9.796e-03 3.633e-03 1.000e+01 0.000e+00
81 133 5.358e-04 2.720e-05 1.168e-02 9.197e-03 3.545e-03 1.000e+01 0.000e+00
82 134 5.062e-04 2.958e-05 1.222e-02 8.745e-03 3.452e-03 1.000e+01 0.000e+00
83 135 4.780e-04 2.820e-05 1.265e-02 8.416e-03 3.365e-03 1.000e+01 0.000e+00
84 136 4.526e-04 2.539e-05 1.295e-02 8.117e-03 3.287e-03 1.000e+01 0.000e+00
85 138 3.653e-04 8.724e-05 3.709e-03 1.832e-02 3.257e-03 1.000e+01 0.000e+00
86 139 4.177e-04 -5.232e-05 1.756e-02 8.603e-03 3.226e-03 1.000e+01 0.000e+00
87 140 3.428e-04 7.488e-05 5.700e-03 4.917e-03 3.193e-03 1.000e+01 3.193e-02
88 141 9.045e-04 -5.617e-04 2.351e-02 1.009e-02 3.066e-03 1.000e+01 3.193e-02
89 142 9.384e-04 -3.393e-05 2.220e-02 1.168e-02 3.083e-03 1.000e+01 3.193e-02
90 143 7.927e-04 1.458e-04 6.223e-03 2.222e-03 3.093e-03 1.000e+01 3.193e-02
91 145 7.819e-04 1.078e-05 5.932e-03 2.430e-03 3.104e-03 1.000e+01 3.193e-02
92 147 7.750e-04 6.824e-06 5.861e-03 2.113e-03 3.105e-03 1.000e+01 3.193e-02
93 149 7.697e-04 5.366e-06 5.770e-03 1.815e-03 3.081e-03 1.000e+01 3.193e-02
94 151 7.670e-04 2.705e-06 5.707e-03 1.488e-03 3.032e-03 1.000e+01 3.193e-02
95 153 7.675e-04 -5.416e-07 5.708e-03 1.365e-03 2.961e-03 1.000e+01 3.193e-02
96 154 1.108e-03 -3.404e-04 2.238e-02 1.492e-02 2.796e-03 1.000e+01 3.193e-02
97 155 7.934e-04 3.145e-04 6.773e-03 1.407e-03 2.583e-03 1.000e+01 3.193e-02
98 156 1.041e-03 -2.475e-04 2.430e-02 1.616e-02 2.354e-03 1.000e+01 3.193e-02
99 157 8.224e-04 2.185e-04 8.995e-03 1.340e-03 2.261e-03 1.000e+01 3.193e-02
100 158 9.843e-04 -1.619e-04 3.040e-02 1.308e-02 2.081e-03 1.000e+01 3.193e-02
101 159 8.262e-04 1.581e-04 8.189e-03 9.091e-04 2.076e-03 1.000e+01 3.193e-02
102 160 8.359e-04 -9.729e-06 3.092e-02 6.702e-03 2.086e-03 1.000e+01 3.193e-02
103 161 8.426e-04 -6.668e-06 2.987e-02 7.654e-03 2.045e-03 1.000e+01 3.193e-02
104 163 8.697e-04 -2.715e-05 2.957e-02 8.741e-03 2.009e-03 1.000e+01 3.193e-02
105 164 8.872e-04 -1.750e-05 2.917e-02 4.917e-02 7.746e-03 1.000e+01 3.193e-02
106 165 8.339e-04 5.337e-05 1.093e-02 8.208e-04 1.950e-03 1.000e+01 3.193e-02
107 166 8.913e-04 -5.740e-05 3.800e-02 1.035e-02 1.902e-03 1.000e+01 3.193e-02
108 167 9.079e-04 -1.664e-05 3.139e-02 1.239e-02 1.820e-03 1.000e+01 3.193e-02
109 168 9.101e-04 -2.214e-06 4.082e-02 1.250e-02 1.827e-03 1.000e+01 3.193e-02
110 169 9.109e-04 -8.115e-07 3.762e-02 1.320e-02 1.850e-03 1.000e+01 3.193e-02
111 170 9.128e-04 -1.916e-06 3.728e-02 1.359e-02 1.875e-03 1.000e+01 3.193e-02
112 171 9.101e-04 2.783e-06 3.648e-02 1.371e-02 1.904e-03 1.000e+01 3.193e-02
113 172 9.041e-04 5.941e-06 3.583e-02 1.355e-02 1.952e-03 1.000e+01 3.193e-02
114 173 8.967e-04 7.398e-06 3.506e-02 1.292e-02 2.050e-03 1.000e+01 3.193e-02
115 174 8.866e-04 1.007e-05 3.422e-02 1.176e-02 2.135e-03 1.000e+01 3.193e-02
116 175 8.736e-04 1.309e-05 3.333e-02 1.014e-02 2.204e-03 1.000e+01 3.193e-02
117 176 8.587e-04 1.486e-05 3.238e-02 1.032e-02 2.260e-03 1.000e+01 3.193e-02
118 177 8.423e-04 1.642e-05 3.135e-02 1.017e-02 2.303e-03 1.000e+01 3.193e-02
119 178 8.221e-04 2.019e-05 3.026e-02 9.528e-03 2.333e-03 1.000e+01 3.193e-02
120 179 7.972e-04 2.492e-05 2.913e-02 8.253e-03 2.368e-03 1.000e+01 3.193e-02
121 180 7.706e-04 2.657e-05 2.797e-02 6.256e-03 2.389e-03 1.000e+01 3.193e-02
122 181 7.496e-04 2.100e-05 2.680e-02 5.684e-03 2.381e-03 1.000e+01 3.193e-02
123 182 6.336e-04 1.160e-04 8.900e-03 7.968e-04 2.415e-03 1.000e+01 3.193e-02
124 184 6.269e-04 6.707e-06 9.023e-03 7.432e-04 2.395e-03 1.000e+01 3.193e-02
125 186 6.207e-04 6.195e-06 8.396e-03 7.892e-04 2.376e-03 1.000e+01 3.193e-02
126 188 6.149e-04 5.764e-06 8.833e-03 7.509e-04 2.372e-03 1.000e+01 3.193e-02
127 190 6.093e-04 5.630e-06 9.512e-03 6.380e-04 2.373e-03 1.000e+01 3.193e-02
128 191 6.660e-04 -5.668e-05 3.988e-02 6.939e-03 2.313e-03 1.000e+01 3.193e-02
129 192 5.866e-04 7.937e-05 7.057e-03 9.666e-04 2.384e-03 1.000e+01 3.193e-02
130 194 5.800e-04 6.620e-06 1.335e-02 6.131e-04 2.388e-03 1.000e+01 3.193e-02
131 196 5.748e-04 5.191e-06 1.232e-02 5.183e-04 2.396e-03 1.000e+01 3.193e-02
132 198 5.698e-04 4.997e-06 1.206e-02 4.834e-04 2.402e-03 1.000e+01 3.193e-02
133 200 5.659e-04 3.937e-06 1.003e-02 3.722e-04 2.410e-03 1.000e+01 3.193e-02
134 201 5.809e-04 -1.505e-05 1.638e-02 4.518e-03 2.442e-03 1.000e+01 3.193e-02
135 202 5.573e-04 2.368e-05 2.406e-03 5.553e-04 2.460e-03 1.000e+01 3.193e-02
136 204 5.559e-04 1.332e-06 2.536e-03 5.522e-04 2.475e-03 1.000e+01 3.193e-02
137 206 5.541e-04 1.855e-06 2.553e-03 5.473e-04 2.488e-03 1.000e+01 3.193e-02
138 208 5.522e-04 1.846e-06 2.565e-03 5.222e-04 2.501e-03 1.000e+01 3.193e-02
139 210 5.504e-04 1.836e-06 2.566e-03 4.946e-04 2.514e-03 1.000e+01 3.193e-02
140 212 5.486e-04 1.794e-06 2.541e-03 4.582e-04 2.526e-03 1.000e+01 3.193e-02
141 214 5.468e-04 1.747e-06 2.542e-03 4.349e-04 2.538e-03 1.000e+01 3.193e-02
142 216 5.451e-04 1.721e-06 2.544e-03 4.190e-04 2.550e-03 1.000e+01 3.193e-02
143 218 5.434e-04 1.699e-06 2.546e-03 4.065e-04 2.561e-03 1.000e+01 3.193e-02
144 220 5.417e-04 1.683e-06 2.547e-03 3.973e-04 2.572e-03 1.000e+01 3.193e-02
145 222 5.401e-04 1.672e-06 2.549e-03 3.909e-04 2.583e-03 1.000e+01 3.193e-02
146 224 5.384e-04 1.666e-06 2.551e-03 3.870e-04 2.594e-03 1.000e+01 3.193e-02
147 226 5.367e-04 1.664e-06 2.553e-03 3.851e-04 2.604e-03 1.000e+01 3.193e-02
148 228 5.351e-04 1.664e-06 2.554e-03 3.851e-04 2.614e-03 1.000e+01 3.193e-02
149 230 5.334e-04 1.667e-06 2.556e-03 3.865e-04 2.624e-03 1.000e+01 3.193e-02
150 232 5.317e-04 1.672e-06 2.557e-03 3.890e-04 2.634e-03 1.000e+01 3.193e-02
151 234 5.301e-04 1.677e-06 2.557e-03 3.925e-04 2.643e-03 1.000e+01 3.193e-02
152 236 5.284e-04 1.682e-06 2.558e-03 3.968e-04 2.652e-03 1.000e+01 3.193e-02
153 238 5.267e-04 1.687e-06 2.558e-03 4.018e-04 2.661e-03 1.000e+01 3.193e-02
154 240 5.250e-04 1.692e-06 2.558e-03 4.073e-04 2.670e-03 1.000e+01 3.193e-02
155 242 5.233e-04 1.696e-06 2.557e-03 4.131e-04 2.678e-03 1.000e+01 3.193e-02
156 244 5.216e-04 1.699e-06 2.557e-03 4.192e-04 2.686e-03 1.000e+01 3.193e-02
157 246 5.199e-04 1.702e-06 2.555e-03 4.255e-04 2.694e-03 1.000e+01 3.193e-02
158 248 5.182e-04 1.704e-06 2.554e-03 4.319e-04 2.701e-03 1.000e+01 3.193e-02
159 250 5.165e-04 1.705e-06 2.552e-03 4.383e-04 2.708e-03 1.000e+01 3.193e-02
160 252 5.148e-04 1.706e-06 2.549e-03 4.446e-04 2.715e-03 1.000e+01 3.193e-02
161 254 5.131e-04 1.706e-06 2.547e-03 4.509e-04 2.721e-03 1.000e+01 3.193e-02
162 256 5.114e-04 1.706e-06 2.544e-03 4.569e-04 2.727e-03 1.000e+01 3.193e-02
163 258 5.097e-04 1.705e-06 2.541e-03 4.628e-04 2.733e-03 1.000e+01 3.193e-02
164 260 5.080e-04 1.703e-06 2.537e-03 4.684e-04 2.738e-03 1.000e+01 3.193e-02
165 262 5.063e-04 1.702e-06 2.533e-03 4.737e-04 2.743e-03 1.000e+01 3.193e-02
166 264 5.046e-04 1.699e-06 2.529e-03 4.787e-04 2.749e-03 1.000e+01 3.193e-02
167 266 5.029e-04 1.697e-06 2.525e-03 4.833e-04 2.763e-03 1.000e+01 3.193e-02
168 268 5.012e-04 1.694e-06 2.522e-03 4.875e-04 2.778e-03 1.000e+01 3.193e-02
169 270 4.995e-04 1.691e-06 2.518e-03 4.913e-04 2.792e-03 1.000e+01 3.193e-02
170 272 4.978e-04 1.687e-06 2.514e-03 4.946e-04 2.807e-03 1.000e+01 3.193e-02
171 274 4.961e-04 1.683e-06 2.510e-03 4.974e-04 2.821e-03 1.000e+01 3.193e-02
172 276 4.944e-04 1.679e-06 2.507e-03 4.998e-04 2.836e-03 1.000e+01 3.193e-02
173 278 4.928e-04 1.674e-06 2.504e-03 5.017e-04 2.851e-03 1.000e+01 3.193e-02
174 280 4.911e-04 1.669e-06 2.501e-03 5.031e-04 2.866e-03 1.000e+01 3.193e-02
175 282 4.894e-04 1.664e-06 2.499e-03 5.040e-04 2.881e-03 1.000e+01 3.193e-02
176 284 4.878e-04 1.658e-06 2.497e-03 5.044e-04 2.896e-03 1.000e+01 3.193e-02
177 286 4.861e-04 1.651e-06 2.496e-03 5.043e-04 2.911e-03 1.000e+01 3.193e-02
178 288 4.845e-04 1.644e-06 2.495e-03 5.037e-04 2.926e-03 1.000e+01 3.193e-02
179 290 4.828e-04 1.637e-06 2.495e-03 5.027e-04 2.941e-03 1.000e+01 3.193e-02
180 292 4.812e-04 1.629e-06 2.495e-03 5.011e-04 2.957e-03 1.000e+01 3.193e-02
181 294 4.796e-04 1.621e-06 2.496e-03 4.991e-04 2.972e-03 1.000e+01 3.193e-02
182 296 4.780e-04 1.612e-06 2.498e-03 4.967e-04 2.987e-03 1.000e+01 3.193e-02
183 298 4.764e-04 1.604e-06 2.499e-03 4.938e-04 3.003e-03 1.000e+01 3.193e-02
184 300 4.748e-04 1.594e-06 2.502e-03 4.905e-04 3.018e-03 1.000e+01 3.193e-02
185 302 4.732e-04 1.585e-06 2.505e-03 4.868e-04 3.033e-03 1.000e+01 3.193e-02
186 304 4.716e-04 1.575e-06 2.508e-03 4.828e-04 3.049e-03 1.000e+01 3.193e-02
187 306 4.701e-04 1.565e-06 2.512e-03 4.785e-04 3.064e-03 1.000e+01 3.193e-02
188 308 4.685e-04 1.554e-06 2.516e-03 4.739e-04 3.079e-03 1.000e+01 3.193e-02
189 310 4.670e-04 1.544e-06 2.521e-03 4.690e-04 3.095e-03 1.000e+01 3.193e-02
190 312 4.654e-04 1.533e-06 2.526e-03 4.639e-04 3.110e-03 1.000e+01 3.193e-02
191 314 4.639e-04 1.521e-06 2.532e-03 4.586e-04 3.125e-03 1.000e+01 3.193e-02
192 316 4.624e-04 1.510e-06 2.537e-03 4.530e-04 3.140e-03 1.000e+01 3.193e-02
193 318 4.609e-04 1.498e-06 2.543e-03 4.474e-04 3.155e-03 1.000e+01 3.193e-02
194 320 4.594e-04 1.485e-06 2.550e-03 4.442e-04 3.169e-03 1.000e+01 3.193e-02
195 322 4.579e-04 1.472e-06 2.556e-03 4.441e-04 3.184e-03 1.000e+01 3.193e-02
196 324 4.565e-04 1.460e-06 2.563e-03 4.438e-04 3.198e-03 1.000e+01 3.193e-02
197 326 4.550e-04 1.446e-06 2.569e-03 4.433e-04 3.213e-03 1.000e+01 3.193e-02
198 328 4.536e-04 1.433e-06 2.576e-03 4.426e-04 3.227e-03 1.000e+01 3.193e-02
199 330 4.522e-04 1.419e-06 2.583e-03 4.417e-04 3.241e-03 1.000e+01 3.193e-02
200 332 4.508e-04 1.405e-06 2.590e-03 4.405e-04 3.254e-03 1.000e+01 3.193e-02
Warning: Maximum number of iterations has been exceeded.
Current function value: 4.508e-04
Constraint violation: 3.254e-03
Total delta_x: 3.533e-01
Iterations: 200
Function evaluations: 332
Jacobian evaluations: 201
Timer: Solution time = 1.21 min
Timer: Avg time per step = 363 ms
==============================================================================================================
Start --> End
Total (sum of squares): 2.946e+01 --> 4.508e-04,
Maximum absolute Quasi-symmetry (1,4) two-term error: 6.229e-01 --> 1.241e-02 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.950e-04 --> 1.025e-06 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.501e-01 --> 6.157e-04 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 5.894e-01 --> 1.174e-02 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.845e-04 --> 9.700e-07 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.420e-01 --> 5.826e-04 (normalized)
Maximum absolute Force error: 4.341e+01 --> 3.398e+03 (N)
Minimum absolute Force error: 4.127e-03 --> 1.121e-01 (N)
Average absolute Force error: 3.601e+00 --> 3.037e+02 (N)
Maximum absolute Force error: 4.261e-05 --> 3.336e-03 (normalized)
Minimum absolute Force error: 4.051e-09 --> 1.101e-07 (normalized)
Average absolute Force error: 3.535e-06 --> 2.982e-04 (normalized)
Aspect ratio: 8.000e+00 --> 8.012e+00 (dimensionless)
Elongation: 1.096e+00 --> 3.000e+00 (dimensionless)
Plasma volume: 3.084e-01 --> 3.085e-01 (m^3)
Plasma volume: 0.000e+00 --> 9.476e-05 (normalized error)
Maximum Rotational transform: 2.477e-01 --> 1.102e+00 (dimensionless)
Minimum Rotational transform: 2.477e-01 --> 1.102e+00 (dimensionless)
Average Rotational transform: 2.477e-01 --> 1.102e+00 (dimensionless)
R boundary error: 0.000e+00 --> 0.000e+00 (m)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
==============================================================================================================
As before, results can be improved by running for more iterations. Note the constraint violation may be larger than desired, so it can be helpful to call eq.solve() at the end to decrease the force error without changing the boundary.
[14]:
eqa.solve();
Building objective: force
Precomputing transforms
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed Psi
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed sheet current
Building objective: self_consistency R
Building objective: self_consistency Z
Building objective: lambda gauge
Building objective: axis R self consistency
Building objective: axis Z self consistency
Number of parameters: 120
Number of objectives: 850
Starting optimization
Using method: lsq-exact
Optimization terminated successfully.
`ftol` condition satisfied.
Current function value: 1.712e-05
Total delta_x: 1.514e-01
Iterations: 3
Function evaluations: 4
Jacobian evaluations: 4
==============================================================================================================
Start --> End
Total (sum of squares): 3.319e-04 --> 1.712e-05,
Maximum absolute Force error: 3.398e+03 --> 7.124e+02 (N)
Minimum absolute Force error: 1.121e-01 --> 5.708e-02 (N)
Average absolute Force error: 3.037e+02 --> 7.105e+01 (N)
Maximum absolute Force error: 4.024e-03 --> 8.437e-04 (normalized)
Minimum absolute Force error: 1.328e-07 --> 6.761e-08 (normalized)
Average absolute Force error: 3.597e-04 --> 8.415e-05 (normalized)
R boundary error: 0.000e+00 --> 2.776e-17 (m)
Z boundary error: 0.000e+00 --> 2.866e-17 (m)
Fixed Psi error: 0.000e+00 --> 0.000e+00 (Wb)
Fixed pressure profile error: 0.000e+00 --> 0.000e+00 (Pa)
Fixed current profile error: 0.000e+00 --> 0.000e+00 (A)
Fixed sheet current error: 0.000e+00 --> 0.000e+00 (~)
==============================================================================================================
From the Boozer plot below we see that we are already doing fairly good for QS:
[15]:
plot_boozer_surface(eqa);
As a final comparison, we can look at the boundary shapes obtained by the different methods. We see that the final shapes are fairly similar, with the proximal method giving slightly more elongation and tighter curvature. We could include additional objectives or constraints to try to reduce this if desired.
[16]:
plot_boundaries(
[eq0, eqa, eqfam[-1]], labels=["Initial", "Augmented Lagrangian", "Proximal"]
);
Further example scripts for precise QS optimization can be found in the desc/examples folder