Advanced 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("../../../"))
[2]:
# from desc import set_device
# set_device("gpu")
[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.12.1+206.g91d82b558,using JAX backend, jax version=0.4.31, jaxlib version=0.4.31, dtype=float64
Using device: CPU, with 13.11 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 = 435 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 17.2 ms
Timer: Objective build = 602 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 40.4 ms
Timer: Objective build = 140 ms
Timer: Proximal projection build = 4.15 sec
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 720 ms
Timer: Linear constraint projection build = 944 ms
Number of parameters: 8
Number of objectives: 460
Timer: Initializing the optimization = 5.86 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 9.402e+01 4.166e+03
1 4 9.113e+01 2.890e+00 3.899e-02 1.897e+04
2 5 7.315e+01 1.798e+01 5.385e-03 8.210e+03
3 6 6.463e+01 8.524e+00 1.423e-02 1.054e+03
4 7 5.889e+01 5.736e+00 2.811e-02 6.835e+02
5 8 5.460e+01 4.292e+00 3.628e-02 1.853e+03
6 9 5.089e+01 3.707e+00 1.541e-02 1.284e+03
7 10 4.674e+01 4.156e+00 4.215e-02 6.953e+03
8 12 4.209e+01 4.645e+00 1.459e-02 4.413e+03
9 13 3.745e+01 4.644e+00 2.158e-02 7.288e+02
10 14 3.585e+01 1.599e+00 2.594e-02 1.018e+04
11 15 3.071e+01 5.135e+00 2.664e-02 6.925e+03
12 16 2.825e+01 2.460e+00 2.315e-02 7.845e+03
13 17 2.734e+01 9.151e-01 2.565e-02 4.245e+03
14 18 2.436e+01 2.977e+00 1.394e-02 3.317e+02
15 43 2.436e+01 0.000e+00 0.000e+00 3.317e+02
Warning: A bad approximation caused failure to predict improvement.
Current function value: 2.436e+01
Total delta_x: 1.528e-01
Iterations: 15
Function evaluations: 43
Jacobian evaluations: 15
Timer: Solution time = 3.07 min
Timer: Avg time per step = 11.5 sec
==============================================================================================================
Start --> End
Total (sum of squares): 9.402e+01 --> 2.436e+01,
Maximum absolute Quasi-symmetry (1,4) two-term error: 6.229e-01 --> 4.064e-01 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.950e-04 --> 9.457e-05 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.501e-01 --> 7.453e-02 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.053e+00 --> 6.871e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 3.296e-04 --> 1.599e-04 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 2.537e-01 --> 1.260e-01 (normalized)
Aspect ratio: 8.000e+00 --> 7.998e+00 (dimensionless)
Maximum absolute Force error: 4.341e+01 --> 1.040e+03 (N)
Minimum absolute Force error: 4.127e-03 --> 8.588e-02 (N)
Average absolute Force error: 3.601e+00 --> 7.963e+01 (N)
Maximum absolute Force error: 6.274e-05 --> 1.503e-03 (normalized)
Minimum absolute Force error: 5.965e-09 --> 1.241e-07 (normalized)
Average absolute Force error: 5.205e-06 --> 1.151e-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)
==============================================================================================================
[8]:
eq2 = run_qh_step(2, eq1)
eqfam.append(eq2)
plot_boozer_surface(eq2);
Building objective: QS two-term
Precomputing transforms
Timer: Precomputing transforms = 37.6 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 17.9 ms
Timer: Objective build = 170 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 38.4 ms
Timer: Objective build = 140 ms
Timer: Proximal projection build = 1.82 sec
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 709 ms
Timer: Linear constraint projection build = 798 ms
Number of parameters: 24
Number of objectives: 460
Timer: Initializing the optimization = 3.37 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.569e+01 1.509e+03
1 5 3.153e+01 4.160e+00 1.565e-02 2.192e+03
2 6 2.921e+01 2.326e+00 1.394e-02 6.361e+03
3 7 2.385e+01 5.353e+00 1.061e-02 5.051e+03
4 8 2.327e+01 5.825e-01 1.961e-02 1.081e+04
5 9 1.729e+01 5.982e+00 4.699e-03 5.886e+02
6 10 1.529e+01 2.001e+00 9.297e-03 1.621e+03
7 11 1.271e+01 2.580e+00 1.943e-02 1.964e+03
8 13 1.116e+01 1.545e+00 1.134e-02 7.326e+02
9 14 9.106e+00 2.057e+00 1.679e-02 1.793e+03
10 15 6.733e+00 2.372e+00 1.743e-02 2.617e+03
11 16 4.927e+00 1.806e+00 2.197e-02 4.138e+03
12 17 2.667e+00 2.260e+00 2.312e-02 2.901e+03
13 19 1.646e+00 1.021e+00 1.257e-02 1.374e+03
14 21 1.340e+00 3.059e-01 6.390e-03 1.872e+02
15 22 1.153e+00 1.872e-01 1.045e-02 1.098e+03
16 24 9.703e-01 1.830e-01 5.438e-03 2.431e+02
17 25 9.057e-01 6.460e-02 1.067e-02 1.204e+03
18 26 8.035e-01 1.022e-01 1.241e-02 1.277e+03
19 28 7.013e-01 1.022e-01 4.387e-03 2.469e+02
20 30 6.838e-01 1.745e-02 2.156e-03 1.099e+02
Warning: Maximum number of iterations has been exceeded.
Current function value: 6.838e-01
Total delta_x: 1.528e-01
Iterations: 20
Function evaluations: 30
Jacobian evaluations: 21
Timer: Solution time = 2.33 min
Timer: Avg time per step = 6.68 sec
==============================================================================================================
Start --> End
Total (sum of squares): 3.569e+01 --> 6.838e-01,
Maximum absolute Quasi-symmetry (1,4) two-term error: 4.064e-01 --> 1.142e-01 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 9.457e-05 --> 2.743e-05 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 7.453e-02 --> 1.174e-02 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 8.317e-01 --> 2.338e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.935e-04 --> 5.614e-05 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.525e-01 --> 2.402e-02 (normalized)
Aspect ratio: 7.998e+00 --> 8.000e+00 (dimensionless)
Maximum absolute Force error: 1.040e+03 --> 1.623e+03 (N)
Minimum absolute Force error: 8.588e-02 --> 1.583e-01 (N)
Average absolute Force error: 7.963e+01 --> 1.267e+02 (N)
Maximum absolute Force error: 1.589e-03 --> 2.480e-03 (normalized)
Minimum absolute Force error: 1.312e-07 --> 2.419e-07 (normalized)
Average absolute Force error: 1.217e-04 --> 1.936e-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 = 37.4 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 21.4 ms
Timer: Objective build = 169 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 41.1 ms
Timer: Objective build = 145 ms
Timer: Proximal projection build = 1.72 sec
Building objective: lcfs R
Building objective: lcfs Z
Building objective: fixed pressure
Building objective: fixed current
Building objective: fixed Psi
Timer: Objective build = 647 ms
Timer: Linear constraint projection build = 747 ms
Number of parameters: 48
Number of objectives: 460
Timer: Initializing the optimization = 3.16 sec
Starting optimization
Using method: proximal-lsq-exact
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 2.213e-01 9.167e+01
1 4 1.674e-01 5.398e-02 6.905e-03 8.037e+02
2 5 1.287e-01 3.862e-02 5.002e-03 6.318e+02
3 6 1.064e-01 2.236e-02 3.663e-03 1.816e+02
4 7 1.010e-01 5.408e-03 4.572e-03 6.785e+02
5 8 7.893e-02 2.203e-02 1.335e-03 8.098e+01
6 9 7.794e-02 9.972e-04 2.410e-03 2.244e+02
7 10 7.051e-02 7.421e-03 7.700e-04 3.033e+01
8 11 7.034e-02 1.699e-04 1.286e-03 5.552e+01
9 12 6.697e-02 3.373e-03 3.802e-04 1.555e+01
10 13 6.563e-02 1.346e-03 6.214e-04 1.007e+01
11 14 6.420e-02 1.421e-03 6.125e-04 9.409e+00
12 36 6.420e-02 0.000e+00 0.000e+00 9.409e+00
Warning: A bad approximation caused failure to predict improvement.
Current function value: 6.420e-02
Total delta_x: 1.872e-02
Iterations: 12
Function evaluations: 36
Jacobian evaluations: 12
Timer: Solution time = 2.48 min
Timer: Avg time per step = 11.4 sec
==============================================================================================================
Start --> End
Total (sum of squares): 2.213e-01 --> 6.420e-02,
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.142e-01 --> 1.021e-01 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 2.743e-05 --> 4.722e-05 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.174e-02 --> 6.202e-03 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.330e-01 --> 1.189e-01 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 3.193e-05 --> 5.496e-05 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 1.366e-02 --> 7.219e-03 (normalized)
Aspect ratio: 8.000e+00 --> 8.000e+00 (dimensionless)
Maximum absolute Force error: 1.623e+03 --> 1.096e+03 (N)
Minimum absolute Force error: 1.583e-01 --> 1.745e-01 (N)
Average absolute Force error: 1.267e+02 --> 8.424e+01 (N)
Maximum absolute Force error: 1.609e-03 --> 1.086e-03 (normalized)
Minimum absolute Force error: 1.569e-07 --> 1.730e-07 (normalized)
Average absolute Force error: 1.256e-04 --> 8.351e-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 = 36.9 ms
Timer: Objective build = 140 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 = 459 ms
Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 39.1 ms
Building objective: aspect ratio
Precomputing transforms
Timer: Precomputing transforms = 18.0 ms
Building objective: elongation
Precomputing transforms
Timer: Precomputing transforms = 19.8 ms
Building objective: volume
Precomputing transforms
Timer: Precomputing transforms = 18.3 ms
Building objective: rotational transform
Precomputing transforms
Timer: Precomputing transforms = 193 ms
Timer: Objective build = 602 ms
Timer: Linear constraint projection build = 1.27 sec
Timer: Linear constraint projection build = 95.7 ms
Number of parameters: 200
Number of objectives: 459
Number of equality constraints: 852
Number of inequality constraints: 2
Timer: Initializing the optimization = 2.76 sec
Starting optimization
Using method: lsq-auglag
Iteration Total nfev Cost Cost reduction Step norm Optimality Constr viol. Penalty param max(|mltplr|)
0 1 9.402e+01 2.528e+04 8.523e-01 1.000e+01 0.000e+00
1 2 8.906e+01 4.955e+00 2.384e-03 2.425e+04 8.528e-01 1.000e+01 0.000e+00
2 3 7.130e+01 1.776e+01 1.018e-02 2.044e+04 8.548e-01 1.000e+01 0.000e+00
3 4 2.553e+01 4.577e+01 4.748e-02 9.383e+03 8.630e-01 1.000e+01 0.000e+00
4 5 7.371e-01 2.479e+01 1.973e-01 1.094e+03 7.379e-01 1.000e+01 0.000e+00
5 6 3.206e-01 4.165e-01 2.321e-01 6.405e+02 6.710e-01 1.000e+01 0.000e+00
6 8 3.905e-01 -6.991e-02 2.024e-01 9.692e+02 6.028e-01 1.000e+01 0.000e+00
7 9 3.293e-01 6.115e-02 1.870e-01 1.082e+03 4.957e-01 1.000e+01 0.000e+00
8 11 2.277e-01 1.016e-01 1.767e-01 7.609e+02 3.931e-01 1.000e+01 0.000e+00
9 13 1.930e-01 3.470e-02 1.457e-01 9.857e+02 2.637e-01 1.000e+01 0.000e+00
10 15 2.033e-01 -1.029e-02 1.671e-01 1.215e+03 1.639e-01 1.000e+01 0.000e+00
11 16 1.208e-01 8.249e-02 1.158e-01 1.192e+03 9.014e-02 1.000e+01 0.000e+00
12 18 1.191e-01 1.690e-03 8.075e-02 1.063e+03 5.598e-02 1.000e+01 0.000e+00
13 19 4.355e-02 7.558e-02 2.406e-02 7.773e+01 5.158e-02 1.000e+01 0.000e+00
14 21 4.172e-02 1.833e-03 1.965e-02 8.652e+01 4.742e-02 1.000e+01 0.000e+00
15 23 3.961e-02 2.108e-03 1.878e-02 8.345e+01 4.361e-02 1.000e+01 0.000e+00
16 25 3.765e-02 1.966e-03 1.992e-02 8.082e+01 4.031e-02 1.000e+01 0.000e+00
17 27 3.573e-02 1.912e-03 1.915e-02 8.138e+01 3.734e-02 1.000e+01 0.000e+00
18 29 3.390e-02 1.831e-03 1.906e-02 8.290e+01 3.463e-02 1.000e+01 0.000e+00
19 31 3.213e-02 1.771e-03 1.906e-02 8.502e+01 3.213e-02 1.000e+01 0.000e+00
20 33 3.041e-02 1.715e-03 1.909e-02 8.796e+01 2.980e-02 1.000e+01 0.000e+00
21 35 2.875e-02 1.661e-03 1.914e-02 9.162e+01 2.821e-02 1.000e+01 0.000e+00
22 37 2.716e-02 1.597e-03 1.918e-02 9.570e+01 2.685e-02 1.000e+01 0.000e+00
23 39 2.564e-02 1.516e-03 2.083e-02 9.968e+01 2.556e-02 1.000e+01 0.000e+00
24 41 2.423e-02 1.411e-03 2.024e-02 1.033e+02 2.433e-02 1.000e+01 0.000e+00
25 43 2.294e-02 1.289e-03 2.023e-02 1.058e+02 2.318e-02 1.000e+01 0.000e+00
26 45 2.178e-02 1.158e-03 2.014e-02 1.071e+02 2.212e-02 1.000e+01 0.000e+00
27 47 2.075e-02 1.032e-03 1.998e-02 1.078e+02 2.116e-02 1.000e+01 0.000e+00
28 49 1.983e-02 9.237e-04 1.974e-02 1.118e+02 2.030e-02 1.000e+01 0.000e+00
29 51 1.899e-02 8.404e-04 1.942e-02 1.133e+02 1.951e-02 1.000e+01 0.000e+00
30 53 1.821e-02 7.825e-04 1.904e-02 1.124e+02 1.907e-02 1.000e+01 0.000e+00
31 55 1.746e-02 7.458e-04 1.860e-02 1.092e+02 1.891e-02 1.000e+01 0.000e+00
32 57 1.673e-02 7.248e-04 1.814e-02 1.044e+02 1.873e-02 1.000e+01 0.000e+00
33 59 1.602e-02 7.143e-04 1.766e-02 9.859e+01 1.851e-02 1.000e+01 0.000e+00
34 61 1.531e-02 7.107e-04 1.717e-02 9.201e+01 1.827e-02 1.000e+01 0.000e+00
35 63 1.460e-02 7.100e-04 1.671e-02 8.642e+01 1.800e-02 1.000e+01 0.000e+00
36 65 1.389e-02 7.096e-04 1.625e-02 8.263e+01 1.771e-02 1.000e+01 0.000e+00
37 67 1.318e-02 7.081e-04 1.580e-02 7.619e+01 1.741e-02 1.000e+01 0.000e+00
38 69 1.248e-02 7.045e-04 1.536e-02 6.708e+01 1.710e-02 1.000e+01 0.000e+00
39 71 1.178e-02 6.967e-04 1.491e-02 5.562e+01 1.686e-02 1.000e+01 0.000e+00
40 73 1.110e-02 6.825e-04 1.444e-02 4.754e+01 1.683e-02 1.000e+01 0.000e+00
41 75 1.044e-02 6.588e-04 1.395e-02 4.004e+01 1.678e-02 1.000e+01 0.000e+00
42 77 9.816e-03 6.231e-04 1.343e-02 3.412e+01 1.671e-02 1.000e+01 0.000e+00
43 79 9.241e-03 5.750e-04 1.287e-02 4.067e+01 1.662e-02 1.000e+01 0.000e+00
44 81 8.724e-03 5.177e-04 1.230e-02 4.680e+01 1.651e-02 1.000e+01 0.000e+00
45 83 8.266e-03 4.573e-04 1.173e-02 5.204e+01 1.637e-02 1.000e+01 0.000e+00
46 85 7.866e-03 4.006e-04 1.119e-02 5.633e+01 1.621e-02 1.000e+01 0.000e+00
47 87 7.514e-03 3.514e-04 1.069e-02 5.968e+01 1.603e-02 1.000e+01 0.000e+00
48 89 7.204e-03 3.101e-04 1.023e-02 6.196e+01 1.585e-02 1.000e+01 0.000e+00
49 91 6.928e-03 2.764e-04 9.752e-03 6.298e+01 1.566e-02 1.000e+01 0.000e+00
50 93 6.676e-03 2.514e-04 9.250e-03 6.276e+01 1.548e-02 1.000e+01 0.000e+00
51 95 6.440e-03 2.361e-04 8.732e-03 6.178e+01 1.531e-02 1.000e+01 0.000e+00
52 97 6.212e-03 2.286e-04 8.238e-03 6.070e+01 1.513e-02 1.000e+01 0.000e+00
53 99 5.987e-03 2.252e-04 7.798e-03 5.993e+01 1.496e-02 1.000e+01 0.000e+00
54 101 5.765e-03 2.216e-04 7.411e-03 5.966e+01 1.478e-02 1.000e+01 0.000e+00
55 103 5.551e-03 2.144e-04 7.054e-03 6.007e+01 1.460e-02 1.000e+01 0.000e+00
56 105 5.347e-03 2.036e-04 6.712e-03 6.121e+01 1.442e-02 1.000e+01 0.000e+00
57 107 5.155e-03 1.916e-04 6.611e-03 6.275e+01 1.424e-02 1.000e+01 0.000e+00
58 109 4.976e-03 1.798e-04 6.294e-03 6.429e+01 1.405e-02 1.000e+01 0.000e+00
59 111 4.807e-03 1.690e-04 6.060e-03 6.548e+01 1.385e-02 1.000e+01 0.000e+00
60 113 4.648e-03 1.588e-04 5.850e-03 6.622e+01 1.366e-02 1.000e+01 0.000e+00
61 115 4.499e-03 1.491e-04 5.665e-03 6.653e+01 1.345e-02 1.000e+01 0.000e+00
62 117 4.359e-03 1.400e-04 5.506e-03 6.647e+01 1.325e-02 1.000e+01 0.000e+00
63 119 4.227e-03 1.315e-04 5.374e-03 6.611e+01 1.304e-02 1.000e+01 0.000e+00
64 121 4.103e-03 1.238e-04 5.269e-03 6.551e+01 1.283e-02 1.000e+01 0.000e+00
65 123 3.986e-03 1.169e-04 5.190e-03 6.473e+01 1.263e-02 1.000e+01 0.000e+00
66 125 3.876e-03 1.108e-04 5.137e-03 6.382e+01 1.242e-02 1.000e+01 0.000e+00
67 127 3.770e-03 1.054e-04 5.105e-03 6.286e+01 1.222e-02 1.000e+01 0.000e+00
68 129 3.670e-03 1.005e-04 5.090e-03 6.191e+01 1.203e-02 1.000e+01 0.000e+00
69 131 3.574e-03 9.607e-05 5.090e-03 6.104e+01 1.184e-02 1.000e+01 0.000e+00
70 133 3.482e-03 9.203e-05 5.102e-03 6.035e+01 1.165e-02 1.000e+01 0.000e+00
71 135 3.393e-03 8.835e-05 5.125e-03 5.992e+01 1.147e-02 1.000e+01 0.000e+00
72 137 3.308e-03 8.495e-05 5.157e-03 5.983e+01 1.129e-02 1.000e+01 0.000e+00
73 139 3.227e-03 8.174e-05 5.198e-03 6.016e+01 1.112e-02 1.000e+01 0.000e+00
74 141 3.148e-03 7.859e-05 5.249e-03 6.097e+01 1.096e-02 1.000e+01 0.000e+00
75 143 3.073e-03 7.539e-05 5.307e-03 6.226e+01 1.079e-02 1.000e+01 0.000e+00
76 145 3.001e-03 7.207e-05 5.370e-03 6.397e+01 1.064e-02 1.000e+01 0.000e+00
77 147 2.932e-03 6.869e-05 5.436e-03 6.598e+01 1.048e-02 1.000e+01 0.000e+00
78 149 2.866e-03 6.541e-05 5.497e-03 6.806e+01 1.033e-02 1.000e+01 0.000e+00
79 151 2.804e-03 6.252e-05 5.547e-03 6.990e+01 1.018e-02 1.000e+01 0.000e+00
80 153 2.744e-03 6.041e-05 5.584e-03 7.118e+01 1.006e-02 1.000e+01 0.000e+00
81 155 2.684e-03 5.945e-05 5.605e-03 7.156e+01 9.930e-03 1.000e+01 0.000e+00
82 157 2.624e-03 5.965e-05 5.615e-03 7.085e+01 9.806e-03 1.000e+01 0.000e+00
83 159 2.564e-03 6.044e-05 5.619e-03 6.898e+01 9.684e-03 1.000e+01 0.000e+00
84 161 2.503e-03 6.095e-05 5.621e-03 6.605e+01 9.564e-03 1.000e+01 0.000e+00
85 163 2.443e-03 6.048e-05 5.623e-03 6.228e+01 9.447e-03 1.000e+01 0.000e+00
86 165 2.384e-03 5.885e-05 5.625e-03 5.789e+01 9.333e-03 1.000e+01 0.000e+00
87 167 2.327e-03 5.624e-05 5.627e-03 5.311e+01 9.219e-03 1.000e+01 0.000e+00
88 169 2.274e-03 5.299e-05 5.628e-03 4.812e+01 9.106e-03 1.000e+01 0.000e+00
89 171 2.225e-03 4.940e-05 5.625e-03 4.308e+01 8.993e-03 1.000e+01 0.000e+00
90 173 2.179e-03 4.574e-05 5.620e-03 3.811e+01 8.878e-03 1.000e+01 0.000e+00
91 175 2.137e-03 4.220e-05 5.611e-03 3.332e+01 8.763e-03 1.000e+01 0.000e+00
92 177 2.098e-03 3.894e-05 5.601e-03 2.878e+01 8.646e-03 1.000e+01 0.000e+00
93 179 2.062e-03 3.606e-05 5.590e-03 2.517e+01 8.527e-03 1.000e+01 0.000e+00
94 181 2.028e-03 3.364e-05 5.579e-03 2.320e+01 8.407e-03 1.000e+01 0.000e+00
95 183 1.997e-03 3.167e-05 5.567e-03 2.246e+01 8.285e-03 1.000e+01 0.000e+00
96 185 1.967e-03 3.015e-05 5.554e-03 2.184e+01 8.163e-03 1.000e+01 0.000e+00
97 187 1.938e-03 2.906e-05 5.541e-03 2.129e+01 8.040e-03 1.000e+01 0.000e+00
98 189 1.909e-03 2.838e-05 5.527e-03 2.116e+01 7.917e-03 1.000e+01 0.000e+00
99 191 1.881e-03 2.807e-05 5.510e-03 2.157e+01 7.794e-03 1.000e+01 0.000e+00
100 193 1.853e-03 2.810e-05 5.490e-03 2.189e+01 7.673e-03 1.000e+01 0.000e+00
101 195 1.825e-03 2.845e-05 5.468e-03 2.207e+01 7.552e-03 1.000e+01 0.000e+00
102 197 1.795e-03 2.909e-05 5.441e-03 2.206e+01 7.434e-03 1.000e+01 0.000e+00
103 199 1.765e-03 2.998e-05 5.409e-03 2.182e+01 7.319e-03 1.000e+01 0.000e+00
104 201 1.734e-03 3.108e-05 5.369e-03 2.133e+01 7.207e-03 1.000e+01 0.000e+00
105 203 1.702e-03 3.228e-05 5.317e-03 2.060e+01 7.099e-03 1.000e+01 0.000e+00
106 205 1.669e-03 3.344e-05 5.249e-03 1.967e+01 6.997e-03 1.000e+01 0.000e+00
107 207 1.634e-03 3.435e-05 5.167e-03 1.943e+01 6.901e-03 1.000e+01 0.000e+00
108 209 1.600e-03 3.475e-05 5.078e-03 1.963e+01 6.813e-03 1.000e+01 0.000e+00
109 211 1.565e-03 3.443e-05 4.997e-03 2.011e+01 6.733e-03 1.000e+01 0.000e+00
110 213 1.532e-03 3.323e-05 4.943e-03 2.119e+01 6.662e-03 1.000e+01 0.000e+00
111 215 1.501e-03 3.120e-05 4.929e-03 2.307e+01 6.599e-03 1.000e+01 0.000e+00
112 217 1.472e-03 2.897e-05 4.960e-03 2.564e+01 6.542e-03 1.000e+01 0.000e+00
113 219 1.443e-03 2.852e-05 5.025e-03 2.839e+01 6.488e-03 1.000e+01 0.000e+00
114 221 1.411e-03 3.270e-05 5.093e-03 3.059e+01 6.437e-03 1.000e+01 0.000e+00
115 223 1.369e-03 4.148e-05 5.125e-03 3.159e+01 6.384e-03 1.000e+01 0.000e+00
116 225 1.320e-03 4.877e-05 5.116e-03 3.067e+01 6.329e-03 1.000e+01 0.000e+00
117 227 1.271e-03 4.882e-05 5.100e-03 2.723e+01 6.270e-03 1.000e+01 0.000e+00
118 229 1.228e-03 4.376e-05 5.118e-03 2.818e+01 6.206e-03 1.000e+01 0.000e+00
119 231 1.189e-03 3.843e-05 5.181e-03 2.777e+01 6.137e-03 1.000e+01 0.000e+00
120 233 1.154e-03 3.495e-05 5.284e-03 2.658e+01 6.065e-03 1.000e+01 0.000e+00
121 235 1.121e-03 3.334e-05 5.516e-03 2.663e+01 5.991e-03 1.000e+01 0.000e+00
122 237 1.088e-03 3.309e-05 5.616e-03 2.722e+01 5.916e-03 1.000e+01 0.000e+00
123 239 1.055e-03 3.330e-05 5.804e-03 2.695e+01 5.842e-03 1.000e+01 0.000e+00
124 241 1.021e-03 3.360e-05 5.989e-03 2.595e+01 5.768e-03 1.000e+01 0.000e+00
125 243 9.875e-04 3.346e-05 6.161e-03 2.442e+01 5.695e-03 1.000e+01 0.000e+00
126 245 9.548e-04 3.273e-05 6.318e-03 2.256e+01 5.623e-03 1.000e+01 0.000e+00
127 247 9.234e-04 3.144e-05 6.456e-03 2.058e+01 5.552e-03 1.000e+01 0.000e+00
128 249 8.937e-04 2.972e-05 6.574e-03 1.862e+01 5.480e-03 1.000e+01 0.000e+00
129 251 8.659e-04 2.773e-05 6.674e-03 1.841e+01 5.407e-03 1.000e+01 0.000e+00
130 253 8.404e-04 2.556e-05 6.762e-03 1.983e+01 5.332e-03 1.000e+01 0.000e+00
131 255 8.170e-04 2.334e-05 6.841e-03 2.102e+01 5.253e-03 1.000e+01 0.000e+00
132 257 7.958e-04 2.124e-05 6.916e-03 2.207e+01 5.172e-03 1.000e+01 0.000e+00
133 259 7.763e-04 1.949e-05 6.987e-03 2.317e+01 5.086e-03 1.000e+01 0.000e+00
134 261 7.580e-04 1.832e-05 7.050e-03 2.456e+01 5.009e-03 1.000e+01 0.000e+00
135 263 7.401e-04 1.791e-05 7.103e-03 2.668e+01 4.977e-03 1.000e+01 0.000e+00
136 265 7.217e-04 1.837e-05 7.140e-03 3.031e+01 4.942e-03 1.000e+01 0.000e+00
137 267 7.022e-04 1.954e-05 1.031e-02 3.706e+01 4.907e-03 1.000e+01 0.000e+00
138 269 6.890e-04 1.317e-05 9.846e-03 5.136e+01 4.872e-03 1.000e+01 0.000e+00
139 271 7.103e-04 -2.129e-05 9.590e-03 8.173e+01 4.853e-03 1.000e+01 0.000e+00
140 272 6.994e-04 1.086e-05 9.091e-03 1.020e+02 4.822e-03 1.000e+01 0.000e+00
141 273 6.545e-04 4.488e-05 8.861e-03 9.806e+01 4.737e-03 1.000e+01 0.000e+00
142 274 6.142e-04 4.031e-05 9.565e-03 8.325e+01 4.647e-03 1.000e+01 0.000e+00
143 275 5.594e-04 5.488e-05 1.044e-02 5.889e+01 4.555e-03 1.000e+01 0.000e+00
144 277 5.275e-04 3.188e-05 1.088e-02 3.439e+01 4.464e-03 1.000e+01 0.000e+00
145 279 5.184e-04 9.089e-06 1.098e-02 2.125e+01 4.372e-03 1.000e+01 0.000e+00
146 281 5.156e-04 2.823e-06 1.090e-02 1.629e+01 4.280e-03 1.000e+01 0.000e+00
147 283 5.137e-04 1.843e-06 1.074e-02 1.508e+01 4.189e-03 1.000e+01 0.000e+00
148 285 5.118e-04 1.894e-06 1.056e-02 1.551e+01 4.099e-03 1.000e+01 0.000e+00
149 286 5.090e-04 2.835e-06 1.520e-02 1.887e+01 4.019e-03 1.000e+01 0.000e+00
150 288 5.039e-04 5.060e-06 1.537e-02 1.428e+01 3.921e-03 1.000e+01 0.000e+00
151 290 5.021e-04 1.789e-06 1.474e-02 1.724e+01 3.834e-03 1.000e+01 0.000e+00
152 292 4.987e-04 3.480e-06 1.462e-02 1.899e+01 3.748e-03 1.000e+01 0.000e+00
153 294 4.943e-04 4.364e-06 1.459e-02 2.074e+01 3.664e-03 1.000e+01 0.000e+00
154 296 4.889e-04 5.417e-06 1.461e-02 2.209e+01 3.586e-03 1.000e+01 0.000e+00
155 298 4.818e-04 7.097e-06 1.795e-02 2.246e+01 3.526e-03 1.000e+01 0.000e+00
156 300 4.742e-04 7.540e-06 2.037e-02 2.347e+01 3.451e-03 1.000e+01 0.000e+00
157 302 4.664e-04 7.815e-06 1.996e-02 2.324e+01 3.381e-03 1.000e+01 0.000e+00
158 304 4.571e-04 9.369e-06 1.947e-02 2.159e+01 3.317e-03 1.000e+01 0.000e+00
159 306 4.494e-04 7.645e-06 2.376e-02 2.276e+01 3.245e-03 1.000e+01 0.000e+00
160 308 4.399e-04 9.530e-06 2.320e-02 2.134e+01 3.174e-03 1.000e+01 0.000e+00
161 310 4.307e-04 9.153e-06 2.336e-02 2.022e+01 3.105e-03 1.000e+01 0.000e+00
162 312 4.216e-04 9.177e-06 2.360e-02 1.892e+01 3.035e-03 1.000e+01 0.000e+00
163 314 4.122e-04 9.332e-06 2.388e-02 1.745e+01 2.966e-03 1.000e+01 0.000e+00
164 316 4.029e-04 9.342e-06 2.415e-02 1.583e+01 2.898e-03 1.000e+01 0.000e+00
165 318 3.936e-04 9.324e-06 2.442e-02 1.436e+01 2.831e-03 1.000e+01 0.000e+00
166 320 3.843e-04 9.274e-06 2.467e-02 1.306e+01 2.765e-03 1.000e+01 0.000e+00
167 322 3.750e-04 9.242e-06 2.491e-02 1.203e+01 2.701e-03 1.000e+01 0.000e+00
168 324 3.658e-04 9.250e-06 2.514e-02 1.088e+01 2.640e-03 1.000e+01 0.000e+00
169 326 3.565e-04 9.296e-06 2.534e-02 9.834e+00 2.581e-03 1.000e+01 0.000e+00
170 328 3.471e-04 9.350e-06 2.552e-02 1.099e+01 2.525e-03 1.000e+01 0.000e+00
171 330 3.378e-04 9.374e-06 2.567e-02 1.175e+01 2.471e-03 1.000e+01 0.000e+00
172 332 3.284e-04 9.336e-06 2.579e-02 1.217e+01 2.420e-03 1.000e+01 0.000e+00
173 334 3.192e-04 9.225e-06 2.587e-02 1.233e+01 2.371e-03 1.000e+01 0.000e+00
174 336 3.102e-04 9.046e-06 2.592e-02 1.231e+01 2.325e-03 1.000e+01 0.000e+00
175 338 3.014e-04 8.813e-06 2.594e-02 1.217e+01 2.281e-03 1.000e+01 0.000e+00
176 340 2.928e-04 8.541e-06 2.592e-02 1.196e+01 2.239e-03 1.000e+01 0.000e+00
177 342 2.846e-04 8.240e-06 2.583e-02 1.215e+01 2.199e-03 1.000e+01 0.000e+00
178 344 2.767e-04 7.910e-06 2.566e-02 1.275e+01 2.161e-03 1.000e+01 0.000e+00
179 346 2.691e-04 7.530e-06 2.533e-02 1.345e+01 2.125e-03 1.000e+01 0.000e+00
180 348 2.621e-04 7.037e-06 2.474e-02 1.440e+01 2.092e-03 1.000e+01 0.000e+00
181 350 2.558e-04 6.252e-06 2.366e-02 1.582e+01 2.062e-03 1.000e+01 0.000e+00
182 352 2.512e-04 4.682e-06 2.167e-02 1.813e+01 2.036e-03 1.000e+01 0.000e+00
183 354 2.500e-04 1.117e-06 1.803e-02 2.186e+01 2.016e-03 1.000e+01 0.000e+00
184 355 2.559e-04 -5.837e-06 1.226e-02 2.945e+01 2.016e-03 1.000e+01 0.000e+00
185 356 2.204e-04 3.544e-05 3.893e-03 3.050e+00 2.017e-03 1.000e+01 0.000e+00
186 358 2.187e-04 1.752e-06 2.642e-03 2.151e+00 2.018e-03 1.000e+01 0.000e+00
187 360 2.174e-04 1.266e-06 1.052e-03 2.396e+00 2.023e-03 1.000e+01 0.000e+00
188 362 2.156e-04 1.800e-06 1.172e-03 2.565e+00 2.029e-03 1.000e+01 0.000e+00
189 364 2.136e-04 2.052e-06 1.264e-03 2.937e+00 2.035e-03 1.000e+01 0.000e+00
190 366 2.114e-04 2.158e-06 1.258e-03 3.021e+00 2.040e-03 1.000e+01 0.000e+00
191 368 2.093e-04 2.081e-06 1.245e-03 3.026e+00 2.044e-03 1.000e+01 0.000e+00
192 370 2.074e-04 1.922e-06 1.293e-03 3.523e+00 2.047e-03 1.000e+01 0.000e+00
193 372 2.057e-04 1.702e-06 1.333e-03 4.030e+00 2.050e-03 1.000e+01 0.000e+00
194 374 2.043e-04 1.446e-06 1.369e-03 4.414e+00 2.052e-03 1.000e+01 0.000e+00
195 376 2.031e-04 1.204e-06 1.404e-03 4.708e+00 2.053e-03 1.000e+01 0.000e+00
196 378 2.021e-04 9.952e-07 1.439e-03 4.952e+00 2.053e-03 1.000e+01 0.000e+00
197 380 2.012e-04 8.186e-07 1.475e-03 5.175e+00 2.053e-03 1.000e+01 0.000e+00
198 382 2.006e-04 6.705e-07 1.511e-03 5.390e+00 2.053e-03 1.000e+01 0.000e+00
199 384 2.000e-04 5.465e-07 1.546e-03 5.599e+00 2.052e-03 1.000e+01 0.000e+00
200 386 1.996e-04 4.432e-07 1.579e-03 5.804e+00 2.051e-03 1.000e+01 0.000e+00
Warning: Maximum number of iterations has been exceeded.
Current function value: 1.996e-04
Constraint violation: 2.051e-03
Total delta_x: 4.197e-01
Iterations: 200
Function evaluations: 386
Jacobian evaluations: 200
Timer: Solution time = 2.61 min
Timer: Avg time per step = 781 ms
==============================================================================================================
Start --> End
Total (sum of squares): 9.402e+01 --> 1.996e-04,
Maximum absolute Quasi-symmetry (1,4) two-term error: 6.229e-01 --> 4.083e-03 (T^3)
Minimum absolute Quasi-symmetry (1,4) two-term error: 1.950e-04 --> 3.880e-08 (T^3)
Average absolute Quasi-symmetry (1,4) two-term error: 1.501e-01 --> 2.301e-04 (T^3)
Maximum absolute Quasi-symmetry (1,4) two-term error: 1.053e+00 --> 6.901e-03 (normalized)
Minimum absolute Quasi-symmetry (1,4) two-term error: 3.296e-04 --> 6.559e-08 (normalized)
Average absolute Quasi-symmetry (1,4) two-term error: 2.537e-01 --> 3.889e-04 (normalized)
Maximum absolute Force error: 4.341e+01 --> 2.243e+03 (N)
Minimum absolute Force error: 4.127e-03 --> 1.218e-01 (N)
Average absolute Force error: 3.601e+00 --> 1.716e+02 (N)
Maximum absolute Force error: 6.274e-05 --> 3.242e-03 (normalized)
Minimum absolute Force error: 5.965e-09 --> 1.761e-07 (normalized)
Average absolute Force error: 5.205e-06 --> 2.480e-04 (normalized)
Aspect ratio: 8.000e+00 --> 7.917e+00 (dimensionless)
Elongation: 1.058e+00 --> 3.000e+00 (dimensionless)
Plasma volume: 3.084e-01 --> 3.084e-01 (m^3)
Plasma volume: 0.000e+00 --> 6.568e-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: 3.112e-06
Total delta_x: 7.051e-02
Iterations: 3
Function evaluations: 4
Jacobian evaluations: 4
==============================================================================================================
Start --> End
Total (sum of squares): 5.430e-05 --> 3.112e-06,
Maximum absolute Force error: 2.243e+03 --> 5.160e+02 (N)
Minimum absolute Force error: 1.218e-01 --> 4.213e-02 (N)
Average absolute Force error: 1.716e+02 --> 4.232e+01 (N)
Maximum absolute Force error: 1.871e-03 --> 4.303e-04 (normalized)
Minimum absolute Force error: 1.016e-07 --> 3.513e-08 (normalized)
Average absolute Force error: 1.431e-04 --> 3.529e-05 (normalized)
R boundary error: 0.000e+00 --> 2.439e-19 (m)
Z boundary error: 0.000e+00 --> 0.000e+00 (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