This page was generated from unit-7-optimization/03_Shape_Derivative_Laplace_SemiAuto.ipynb.

7.3 PDE-Constrained Shape Optimization (semi-automated)

We want to solve the PDE-constrained shape optimization problem

\[\underset{\Omega\subset \mathsf{D}}{\mbox{min}} \; J(u) := \int_\Omega |u-u_d|^q \; dx, \quad q\ge 2\]

subject to that \((\Omega,u)\) satisfy

\[\int_\Omega \nabla u \cdot \nabla v \; dx = \int_\Omega f v \; dx \; \quad \text{ for all } v \in H_0^1(\Omega),\]

where \(\Omega \subset \mathbb R^2\) for given \(u_d, \, f \in C^1(\mathbb R^2)\).

Here, we want to compute the shape derivative by differentiation of a suitably defined perturbed Lagrangian using automated differentiation. For details we refer to

16. Gangl, K. Sturm, M. Neunteufel, J. Schöberl. Fully and Semi-Automated Shape Differentiation in NGSolve, Struct. Multidisc. Optim., 63, pp.1579-1607, 2021.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry

geo = SplineGeometry()
geo.AddCircle(c=(0.5,0.5), r=0.5, bc = 'circle')
mesh = Mesh(geo.GenerateMesh(maxh = 0.08))

#given data of our problem (chosen such that \Omega^* = [0,1]^2 is the optimal shape)
f = CoefficientFunction(2*y*(1-y)+2*x*(1-x))
ud = x*(1-x)*y*(1-y)

grad_f = CoefficientFunction( (f.Diff(x), f.Diff(y) ) )
grad_ud = CoefficientFunction( (ud.Diff(x), ud.Diff(y) ) )

State equation

fes = H1(mesh, order=2, dirichlet=".*")
u, v = fes.TnT()
gfu = GridFunction(fes)
scene = Draw(gfu, mesh, "state")

a = BilinearForm(fes, symmetric=True)
a += grad(u)*grad(v)*dx

fstate = LinearForm(fes)
fstate += f*v*dx

def SolveStateEquation():
    rhs = gfu.vec.CreateVector()
    rhs.data = fstate.vec - a.mat * gfu.vec
    update = gfu.vec.CreateVector()
    update.data = a.mat.Inverse(fes.FreeDofs()) * rhs
    gfu.vec.data += update

a.Assemble()
fstate.Assemble()
SolveStateEquation()

scene.Redraw()

Adjoint equation

We set up the adjoint equation

\[\mbox{Find } p \in H_0^1(\Omega): \int_\Omega \nabla w \cdot \nabla p \, \mbox dx = - \partial_u J(u)(w) \quad \text{ for all } w \in H_0^1(\Omega)\]

where \(u\) is the solution to the state equation. For \(J(u) = \int_\Omega |u-u_d|^q \mbox dx\), we get

\[\partial_u J(u)(w) = q \int_\Omega (u-u_d)^{d-1}w \,\mbox dx.\]

However, we can also use the Diff(…) command:

q=4
def Cost(u):
    return (u-ud)**q*dx

p, w = fes.TnT()
gfp = GridFunction(fes)
scene = Draw (gfp, mesh, "adjoint")

fadjoint = LinearForm(fes)
fadjoint += -1*Cost(gfu).Diff(gfu,w)

def SolveAdjointEquation():
    rhs = gfp.vec.CreateVector()
    rhs.data = fadjoint.vec - a.mat.T * gfp.vec
    update = gfp.vec.CreateVector()
    update.data = a.mat.Inverse(fes.FreeDofs()).T * rhs
    gfp.vec.data += update

fadjoint.Assemble()
SolveAdjointEquation()
scene.Redraw()

Note that (for linear problems) the operator on the left hand side of the adjoint equation is just the transpose of the state operator.

Automatic Shape Differentiation

The formula for the shape derivative was derived as the partial derivative of the perturbed Lagrangian (brought back to the original domain):

\[d\mathcal J(\Omega; X) = \frac{\partial}{\partial t} \left. \left( \int_\Omega \xi(t)|u - u_d^t|^q \mbox{d} x + \int_{\Omega} (F_t^{-\top}\nabla u) \cdot (F_t^{-\top} \nabla p) \xi(t) \, dx - \int_{\Omega} \xi(t) f^t p \,dx \right) \right\rvert_{t=0}\]

where

• \(T_t(x)=x+tX(x)=y\)

• \(F_t = DT_t = I+t DX\)

• \(\xi(t) = \mbox{det}(F_t)\)

• \(u_d^t = u_d \circ T_t\)

• \(f^t = f \circ T_t\)

The integrand depends on the parameter \(t\) only via \(\xi(t)\), \(F_t\) and \(T_t\). Denoting the integrand by \(G^{u,p}\), the derivative is given by

\[\begin{split}\begin{array}{rl} d\mathcal J(\Omega; X) =& \frac{d G^{u,p}}{d \xi} \frac{d \xi}{d t} + \frac{d G^{u,p}}{d F} \frac{d F}{d t} + \frac{d G^{u,p}}{dy} \cdot \frac{d T_t}{dt} \\ =& \frac{d G^{u,p}}{d \xi} \mbox{div}(X) + \frac{d G^{u,p}}{d F} DX + \frac{d G^{u,p}}{dy} \cdot X \end{array}\end{split}\]

We define the perturbed PDE and cost function involving parameters \(F\) and \(J\) representing the Jacobian of a shape transformation and its determinant, respectively. While their values are set to unity and thus does not influence the solution of the state or adjoint equation, this procedure allows us to differentiate with respect to them.

J = Parameter(1)
F = Id(2)
#Finv = 2*Id(2) - F    # consistent linearization at F = I (avoids calling Inv() for speedup)

def Equation(u,v):
    return ( (Inv(F.trans)*grad(u))*(Inv(F.trans)*grad(v))-f*v)*J*dx
#    return ( (Finv.trans*grad(u))*(Finv.trans*grad(v))-f*v)*J*dx

def CostAuto(u):
    return J*(u-ud)**q*dx

VEC = H1(mesh, order=2, dim=2)
PHI, X = VEC.TnT()

Lagrangian = CostAuto(gfu) + Equation(gfu,gfp)
dJOmegaAuto = LinearForm(VEC)
dJOmegaAuto += Lagrangian.Diff(J, div(X))
dJOmegaAuto += Lagrangian.Diff(F, grad(X).trans)   # grad(X).trans is Jacobian
dJOmegaAuto += Lagrangian.Diff(x, X[0]) + Lagrangian.Diff(y, X[1])

b = BilinearForm(VEC)
b += InnerProduct(grad(X),grad(PHI))*dx + InnerProduct(X,PHI)*dx

gfX = GridFunction(VEC)

def SolveDeformationEquationAuto():
    rhs = gfX.vec.CreateVector()
    rhs.data = dJOmegaAuto.vec - b.mat * gfX.vec
    update = gfX.vec.CreateVector()
    update.data = b.mat.Inverse(VEC.FreeDofs()) * rhs
    gfX.vec.data += update

b.Assemble()
dJOmegaAuto.Assemble()
SolveDeformationEquationAuto()
Draw(-gfX, mesh, "gfX")

BaseWebGuiScene

# gfset denotes the deformation of the original domain and will be updated during the shape optimization
gfset = GridFunction(VEC)
gfset.Set((0,0))
mesh.SetDeformation(gfset)
sceneSet = Draw(gfset,mesh,"gfset")
SetVisualization (deformation=True)

gfset.Set((0,0))
mesh.SetDeformation(gfset)
print('Cost at initial design', Integrate (CostAuto(gfu), mesh))

scale = 0.5 / Norm(gfX.vec)
gfset.vec.data -= scale * gfX.vec
mesh.SetDeformation(gfset)
sceneSet.Redraw()

Cost at initial design 5.152244130814951e-09

a.Assemble()
fstate.Assemble()
SolveStateEquation()
print('Cost at new design', Integrate (CostAuto(gfu), mesh))

Cost at new design 9.793644493069302e-10

Equation can also be used to define the bilinear form. The following defines the same bilinear form as above:

aAuto = BilinearForm(fes, symmetric=True)
aAuto += Equation(u,v)

Thus, the user has to enter the PDE (in its transformed form) only once.

Finally, let us again run the full algorithm:

#reset to and solve for initial configuration
gfset.Set((0,0))
mesh.SetDeformation(gfset)
scene = Draw(gfset,mesh,"gfset")
SetVisualization (deformation=True)

a.Assemble()
fstate.Assemble()
SolveStateEquation()

LineSearch = False

iter_max = 600
Jold = Integrate(CostAuto(gfu), mesh)
converged = False
for k in range(iter_max):
    scene.Redraw()
    print('cost at iteration', k, ': ', Jold)
    mesh.SetDeformation(gfset)

    a.Assemble()
    fstate.Assemble()
    SolveStateEquation()

    fadjoint.Assemble()
    SolveAdjointEquation()

    b.Assemble()
    dJOmegaAuto.Assemble()
    SolveDeformationEquationAuto()

    scale = 0.01 / Norm(gfX.vec)
    gfsetOld = gfset
    gfset.vec.data -= scale * gfX.vec

    Jnew = Integrate(CostAuto(gfu), mesh)

    if LineSearch:
        while Jnew > Jold and scale > 1e-12:
            scale = scale / 2

            if scale <= 1e-12:
                converged = True
                break

            gfset.vec.data = gfsetOld.vec - scale * gfX.vec
            mesh.SetDeformation(gfset)

            a.Assemble()
            fstate.Assemble()
            SolveStateEquation()
            Jnew = Integrate(CostAuto(gfu), mesh)

    if converged==True:
        break
    Jold = Jnew

cost at iteration 0 :  5.152244130814932e-09
cost at iteration 1 :  4.942040315052465e-09
cost at iteration 2 :  4.5850808690470135e-09
cost at iteration 3 :  4.249439885598259e-09
cost at iteration 4 :  3.934229274496262e-09
cost at iteration 5 :  3.6385740544274783e-09
cost at iteration 6 :  3.3616132927088307e-09
cost at iteration 7 :  3.102501038477676e-09
cost at iteration 8 :  2.8604072470352507e-09
cost at iteration 9 :  2.634518692682722e-09
cost at iteration 10 :  2.424039866929757e-09
cost at iteration 11 :  2.2281938583607805e-09
cost at iteration 12 :  2.0462232096623734e-09
cost at iteration 13 :  1.877390746276767e-09
cost at iteration 14 :  1.7209803697550375e-09
cost at iteration 15 :  1.576297807005677e-09
cost at iteration 16 :  1.4426713040877673e-09
cost at iteration 17 :  1.3194522497306325e-09
cost at iteration 18 :  1.2060157090339114e-09
cost at iteration 19 :  1.1017608413450996e-09
cost at iteration 20 :  1.0061111674975642e-09
cost at iteration 21 :  9.185146395774446e-10
cost at iteration 22 :  8.384434500547998e-10
cost at iteration 23 :  7.653934950250658e-10
cost at iteration 24 :  6.988833766917684e-10
cost at iteration 25 :  6.384527911306818e-10
cost at iteration 26 :  5.836600972575326e-10
cost at iteration 27 :  5.340788023199241e-10
cost at iteration 28 :  4.892926352905075e-10
cost at iteration 29 :  4.488888372771933e-10
cost at iteration 30 :  4.124493475644511e-10
cost at iteration 31 :  3.795398821676597e-10
cost at iteration 32 :  3.496978994973505e-10
cost at iteration 33 :  3.224229646055943e-10
cost at iteration 34 :  2.9717833433569407e-10
cost at iteration 35 :  2.734208497181494e-10
cost at iteration 36 :  2.506798047504566e-10
cost at iteration 37 :  2.2867865248925986e-10
cost at iteration 38 :  2.0741877930875087e-10
cost at iteration 39 :  1.8711564698267465e-10
cost at iteration 40 :  1.680225990721038e-10
cost at iteration 41 :  1.5030409913780633e-10
cost at iteration 42 :  1.340169892860871e-10
cost at iteration 43 :  1.1914272597739713e-10
cost at iteration 44 :  1.056211889134839e-10
cost at iteration 45 :  9.337291061389827e-11
cost at iteration 46 :  8.231145824835373e-11
cost at iteration 47 :  7.23497539531662e-11
cost at iteration 48 :  6.340308482728872e-11
cost at iteration 49 :  5.539042161627136e-11
cost at iteration 50 :  4.823491672945069e-11
cost at iteration 51 :  4.186402390638018e-11
cost at iteration 52 :  3.620945390576394e-11
cost at iteration 53 :  3.12070653802054e-11
cost at iteration 54 :  2.6796734440610483e-11
cost at iteration 55 :  2.2922220688921818e-11
cost at iteration 56 :  1.9531036130474335e-11
cost at iteration 57 :  1.6574318562732776e-11
cost at iteration 58 :  1.4006709060443762e-11
cost at iteration 59 :  1.1786232362324965e-11
cost at iteration 60 :  9.874178616578768e-12
cost at iteration 61 :  8.234984789862004e-12
cost at iteration 62 :  6.8361140194980614e-12
cost at iteration 63 :  5.647931255556368e-12
cost at iteration 64 :  4.643573755009424e-12
cost at iteration 65 :  3.79881527448045e-12
cost at iteration 66 :  3.0919233641335144e-12
cost at iteration 67 :  2.5035095464609645e-12
cost at iteration 68 :  2.016373380151357e-12
cost at iteration 69 :  1.6153408584893742e-12
cost at iteration 70 :  1.2871018492123458e-12
cost at iteration 71 :  1.0200425660121514e-12
cost at iteration 72 :  8.04101811907473e-13
cost at iteration 73 :  6.306625713390689e-13
cost at iteration 74 :  4.937810587536974e-13
cost at iteration 75 :  4.038743115985516e-13
cost at iteration 76 :  3.6977366617695585e-13
cost at iteration 77 :  3.391503989858549e-13
cost at iteration 78 :  3.0994728647082915e-13
cost at iteration 79 :  2.8418809614113074e-13
cost at iteration 80 :  2.602015679563889e-13
cost at iteration 81 :  2.3870400427207476e-13
cost at iteration 82 :  2.189402326498232e-13
cost at iteration 83 :  2.0107235167501253e-13
cost at iteration 84 :  1.847560962827504e-13
cost at iteration 85 :  1.6993450676810084e-13
cost at iteration 86 :  1.5644707213282408e-13
cost at iteration 87 :  1.441612852898547e-13
cost at iteration 88 :  1.3299888316163557e-13
cost at iteration 89 :  1.2281207207771698e-13
cost at iteration 90 :  1.1355919196768976e-13
cost at iteration 91 :  1.0510144282710361e-13
cost at iteration 92 :  9.741322745022404e-14
cost at iteration 93 :  9.037356222322941e-14
cost at iteration 94 :  8.396375795060642e-14
cost at iteration 95 :  7.808287254019635e-14
cost at iteration 96 :  7.271510922055711e-14
cost at iteration 97 :  6.77789668727767e-14
cost at iteration 98 :  6.325978293235702e-14
cost at iteration 99 :  5.909367708276851e-14
cost at iteration 100 :  5.5266312266776365e-14
cost at iteration 101 :  5.1729119032482634e-14
cost at iteration 102 :  4.8467672516927466e-14
cost at iteration 103 :  4.5446296325404274e-14
cost at iteration 104 :  4.265025279546492e-14
cost at iteration 105 :  4.005445437800318e-14
cost at iteration 106 :  3.764374632744396e-14
cost at iteration 107 :  3.540156117394774e-14
cost at iteration 108 :  3.331231820891711e-14
cost at iteration 109 :  3.1366187136757206e-14
cost at iteration 110 :  2.954722865704113e-14
cost at iteration 111 :  2.7850855087319455e-14
cost at iteration 112 :  2.6260891448208297e-14
cost at iteration 113 :  2.4776761800528488e-14
cost at iteration 114 :  2.3382213446773128e-14
cost at iteration 115 :  2.2079678770587107e-14
cost at iteration 116 :  2.0853002877993295e-14
cost at iteration 117 :  1.9706812337217115e-14
cost at iteration 118 :  1.862523000766937e-14
cost at iteration 119 :  1.7614416588232597e-14
cost at iteration 120 :  1.6658948127921077e-14
cost at iteration 121 :  1.5765982826450304e-14
cost at iteration 122 :  1.4920714712236086e-14
cost at iteration 123 :  1.4130859202355867e-14
cost at iteration 124 :  1.3382376270812139e-14
cost at iteration 125 :  1.2683167514497975e-14
cost at iteration 126 :  1.202007901338602e-14
cost at iteration 127 :  1.1400864336218601e-14
cost at iteration 128 :  1.0813330171009113e-14
cost at iteration 129 :  1.0264752199471035e-14
cost at iteration 130 :  9.743925711461055e-15
cost at iteration 131 :  9.257338158356138e-15
cost at iteration 132 :  8.794835075022783e-15
cost at iteration 133 :  8.361938256590832e-15
cost at iteration 134 :  7.94975512261959e-15
cost at iteration 135 :  7.562856884085308e-15
cost at iteration 136 :  7.193878599785462e-15
cost at iteration 137 :  6.846491064372645e-15
cost at iteration 138 :  6.514948805575775e-15
cost at iteration 139 :  6.201993142147536e-15
cost at iteration 140 :  5.903399994243567e-15
cost at iteration 141 :  5.6209360772933205e-15
cost at iteration 142 :  5.351731791472277e-15
cost at iteration 143 :  5.096582419347018e-15
cost at iteration 144 :  4.853808998451093e-15
cost at iteration 145 :  4.623288548177376e-15
cost at iteration 146 :  4.4043868182233105e-15
cost at iteration 147 :  4.1961450156421914e-15
cost at iteration 148 :  3.998842878687953e-15
cost at iteration 149 :  3.8107802788536896e-15
cost at iteration 150 :  3.633031539817425e-15
cost at iteration 151 :  3.4632519214800258e-15
cost at iteration 152 :  3.3031996120242764e-15
cost at iteration 153 :  3.1499799330321967e-15
cost at iteration 154 :  3.0059307879100035e-15
cost at iteration 155 :  2.867699581347287e-15
cost at iteration 156 :  2.738103985676431e-15
cost at iteration 157 :  2.6134257441728256e-15
cost at iteration 158 :  2.4968635516442203e-15
cost at iteration 159 :  2.38443555774674e-15
cost at iteration 160 :  2.279621895111495e-15
cost at iteration 161 :  2.1783221309662862e-15
cost at iteration 162 :  2.0842133225064412e-15
cost at iteration 163 :  1.993409186588248e-15
cost at iteration 164 :  1.909848596847841e-15
cost at iteration 165 :  1.8308844456745757e-15
cost at iteration 166 :  1.7603512878472203e-15
cost at iteration 167 :  1.6953726890043395e-15
cost at iteration 168 :  1.6337580822264638e-15
cost at iteration 169 :  1.5755635686952429e-15
cost at iteration 170 :  1.5186186097079618e-15
cost at iteration 171 :  1.465086425287961e-15
cost at iteration 172 :  1.4125312442929655e-15
cost at iteration 173 :  1.3633234413160113e-15
cost at iteration 174 :  1.3151003514189375e-15
cost at iteration 175 :  1.2700675006280674e-15
cost at iteration 176 :  1.226032639454432e-15
cost at iteration 177 :  1.1849951382265766e-15
cost at iteration 178 :  1.144936271885594e-15
cost at iteration 179 :  1.1076635467262861e-15
cost at iteration 180 :  1.0713194283535159e-15
cost at iteration 181 :  1.0375421921162177e-15
cost at iteration 182 :  1.0046226390029162e-15
cost at iteration 183 :  9.740506483991214e-16
cost at iteration 184 :  9.442543155934123e-16
cost at iteration 185 :  9.165922758378797e-16
cost at iteration 186 :  8.896205658754849e-16
cost at iteration 187 :  8.645802600791525e-16
cost at iteration 188 :  8.401471777233702e-16
cost at iteration 189 :  8.174557391862637e-16
cost at iteration 190 :  7.95294172079007e-16
cost at iteration 191 :  7.746990604426485e-16
cost at iteration 192 :  7.545642017578473e-16
cost at iteration 193 :  7.358356189500002e-16
cost at iteration 194 :  7.175062626986021e-16
cost at iteration 195 :  7.0043772220014e-16
cost at iteration 196 :  6.83716077240544e-16
cost at iteration 197 :  6.681237361936978e-16
cost at iteration 198 :  6.528342847574644e-16
cost at iteration 199 :  6.385555190391363e-16
cost at iteration 200 :  6.245433219296647e-16
cost at iteration 201 :  6.114348986356576e-16
cost at iteration 202 :  5.98563637869376e-16
cost at iteration 203 :  5.864997317207091e-16
cost at iteration 204 :  5.746496889312711e-16
cost at iteration 205 :  5.635199036823379e-16
cost at iteration 206 :  5.525860006247824e-16
cost at iteration 207 :  5.422934930100947e-16
cost at iteration 208 :  5.321834675266953e-16
cost at iteration 209 :  5.226432264394739e-16
cost at iteration 210 :  5.132759802960449e-16
cost at iteration 211 :  5.04413283780636e-16
cost at iteration 212 :  4.957174143768553e-16
cost at iteration 213 :  4.87466467852465e-16
cost at iteration 214 :  4.7937898066731075e-16
cost at iteration 215 :  4.716817283990702e-16
cost at iteration 216 :  4.641469183892862e-16
cost at iteration 217 :  4.569520140990757e-16
cost at iteration 218 :  4.499204999808767e-16
cost at iteration 219 :  4.431824198080625e-16
cost at iteration 220 :  4.3661031365091153e-16
cost at iteration 221 :  4.3028859413004123e-16
cost at iteration 222 :  4.2413678887625054e-16
cost at iteration 223 :  4.181953733069096e-16
cost at iteration 224 :  4.124289319615828e-16
cost at iteration 225 :  4.068356099179253e-16
cost at iteration 226 :  4.0142324174716734e-16
cost at iteration 227 :  3.9614916813659984e-16
cost at iteration 228 :  3.910627789680601e-16
cost at iteration 229 :  3.8608206077185916e-16
cost at iteration 230 :  3.8129636622816025e-16
cost at iteration 231 :  3.765857067084717e-16
cost at iteration 232 :  3.7207789882330614e-16
cost at iteration 233 :  3.6761629049167306e-16
cost at iteration 234 :  3.6336574961092133e-16
cost at iteration 235 :  3.5913420859252025e-16
cost at iteration 236 :  3.5512225373229503e-16
cost at iteration 237 :  3.5110358932165345e-16
cost at iteration 238 :  3.4731326124524073e-16
cost at iteration 239 :  3.434918754410022e-16
cost at iteration 240 :  3.399077476394114e-16
cost at iteration 241 :  3.3626946028307666e-16
cost at iteration 242 :  3.3287747381834904e-16
cost at iteration 243 :  3.2940936966289136e-16
cost at iteration 244 :  3.261966884789644e-16
cost at iteration 245 :  3.2288698309760543e-16
cost at iteration 246 :  3.198418669398781e-16
cost at iteration 247 :  3.166797888701883e-16
cost at iteration 248 :  3.1379148139979635e-16
cost at iteration 249 :  3.1076716832069005e-16
cost at iteration 250 :  3.0802579837855145e-16
cost at iteration 251 :  3.0513020545184464e-16
cost at iteration 252 :  3.025266997344011e-16
cost at iteration 253 :  2.9975151852030033e-16
cost at iteration 254 :  2.972775241844618e-16
cost at iteration 255 :  2.9461511077193494e-16
cost at iteration 256 :  2.9226292670054327e-16
cost at iteration 257 :  2.897062378873735e-16
cost at iteration 258 :  2.8746875352619843e-16
cost at iteration 259 :  2.850112900464336e-16
cost at iteration 260 :  2.828819308751007e-16
cost at iteration 261 :  2.805176868089569e-16
cost at iteration 262 :  2.784903656362595e-16
cost at iteration 263 :  2.7621378325412945e-16
cost at iteration 264 :  2.7428285663742125e-16
cost at iteration 265 :  2.720887860286002e-16
cost at iteration 266 :  2.7024901520951365e-16
cost at iteration 267 :  2.6813267813107392e-16
cost at iteration 268 :  2.663791939596475e-16
cost at iteration 269 :  2.643361514131085e-16
cost at iteration 270 :  2.6266442280042437e-16
cost at iteration 271 :  2.6069054590670973e-16
cost at iteration 272 :  2.590963514047489e-16
cost at iteration 273 :  2.5718779520142113e-16
cost at iteration 274 :  2.5566719735941047e-16
cost at iteration 275 :  2.5382037719069413e-16
cost at iteration 276 :  2.52369699381016e-16
cost at iteration 277 :  2.5058126959174554e-16
cost at iteration 278 :  2.4919707503579335e-16
cost at iteration 279 :  2.474639097152262e-16
cost at iteration 280 :  2.4614298247273943e-16
cost at iteration 281 :  2.4446215802470727e-16
cost at iteration 282 :  2.4320148573822024e-16
cost at iteration 283 :  2.4157026508067735e-16
cost at iteration 284 :  2.403670232913649e-16
cost at iteration 285 :  2.3878284151190023e-16
cost at iteration 286 :  2.3763437938453045e-16
cost at iteration 287 :  2.360948306982518e-16
cost at iteration 288 :  2.3499865801179527e-16
cost at iteration 289 :  2.335014838860535e-16
cost at iteration 290 :  2.324552591626938e-16
cost at iteration 291 :  2.3099833748874366e-16
cost at iteration 292 :  2.299998571481534e-16
cost at iteration 293 :  2.285811923534648e-16
cost at iteration 294 :  2.276283807921332e-16
cost at iteration 295 :  2.2624609479886886e-16
cost at iteration 296 :  2.253369953046784e-16
cost at iteration 297 :  2.239893192508544e-16
cost at iteration 298 :  2.231220856732583e-16
cost at iteration 299 :  2.218073523219687e-16
cost at iteration 300 :  2.2098024142623818e-16
cost at iteration 301 :  2.1969687819681792e-16
cost at iteration 302 :  2.18908242638487e-16
cost at iteration 303 :  2.1765476520047076e-16
cost at iteration 304 :  2.1690304706279428e-16
cost at iteration 305 :  2.1567805344012128e-16
cost at iteration 306 :  2.1496177828295586e-16
cost at iteration 307 :  2.1376394342118567e-16
cost at iteration 308 :  2.1308171479525445e-16
cost at iteration 309 :  2.119097855499005e-16
cost at iteration 310 :  2.1126027993476103e-16
cost at iteration 311 :  2.101130704428174e-16
cost at iteration 312 :  2.0949503257106578e-16
cost at iteration 313 :  2.083714199720127e-16
cost at iteration 314 :  2.0778365850600407e-16
cost at iteration 315 :  2.0668257898182499e-16
cost at iteration 316 :  2.0612396251226557e-16
cost at iteration 317 :  2.0504440761911728e-16
cost at iteration 318 :  2.0451386095812034e-16
cost at iteration 319 :  2.0345487422497733e-16
cost at iteration 320 :  2.029513749684247e-16
cost at iteration 321 :  2.0191204874046022e-16
cost at iteration 322 :  2.0143462407727032e-16
cost at iteration 323 :  2.0041409658377315e-16
cost at iteration 324 :  1.9996182033141956e-16
cost at iteration 325 :  1.9895927296007836e-16
cost at iteration 326 :  1.9853126280792777e-16
cost at iteration 327 :  1.9754591756882086e-16
cost at iteration 328 :  1.9714133251236317e-16
cost at iteration 329 :  1.961724496768217e-16
cost at iteration 330 :  1.957904876273831e-16
cost at iteration 331 :  1.9483736352783183e-16
cost at iteration 332 :  1.9447725908397684e-16
cost at iteration 333 :  1.9353922406242817e-16
cost at iteration 334 :  1.932002464302995e-16
cost at iteration 335 :  1.922766629240335e-16
cost at iteration 336 :  1.9195811397504393e-16
cost at iteration 337 :  1.910483747291865e-16
cost at iteration 338 :  1.9074958718456637e-16
cost at iteration 339 :  1.8985311358198643e-16
cost at iteration 340 :  1.8957344931455586e-16
cost at iteration 341 :  1.886896898143024e-16
cost at iteration 342 :  1.884285382587229e-16
cost at iteration 343 :  1.875569669351599e-16
cost at iteration 344 :  1.8731374359860514e-16
cost at iteration 345 :  1.864538587737831e-16
cost at iteration 346 :  1.862280038396937e-16
cost at iteration 347 :  1.8537932680217822e-16
cost at iteration 348 :  1.8517030382042721e-16
cost at iteration 349 :  1.843323776245457e-16
cost at iteration 350 :  1.8413967228174857e-16
cost at iteration 351 :  1.8331206062132077e-16
cost at iteration 352 :  1.8313517958576101e-16
cost at iteration 353 :  1.823174657372084e-16
cost at iteration 354 :  1.821559355731144e-16
cost at iteration 355 :  1.8134772140295692e-16
cost at iteration 356 :  1.8120108754933917e-16
cost at iteration 357 :  1.8040199258169182e-16
cost at iteration 358 :  1.8026981839143316e-16
cost at iteration 359 :  1.7947947893124026e-16
cost at iteration 360 :  1.793613447663281e-16
cost at iteration 361 :  1.7857941307458443e-16
cost at iteration 362 :  1.7847491545375535e-16
cost at iteration 363 :  1.7770105897118833e-16
cost at iteration 364 :  1.7760980976654561e-16
cost at iteration 365 :  1.7684371038233778e-16
cost at iteration 366 :  1.7676533606184192e-16
cost at iteration 367 :  1.7600668942448574e-16
cost at iteration 368 :  1.759408303372099e-16
cost at iteration 369 :  1.7518934520467335e-16
cost at iteration 370 :  1.7513565490622787e-16
cost at iteration 371 :  1.7439105253269426e-16
cost at iteration 372 :  1.743491971483068e-16
cost at iteration 373 :  1.736112107052476e-16
cost at iteration 374 :  1.7358086832798953e-16
cost at iteration 375 :  1.7284924235719093e-16
cost at iteration 376 :  1.7283010247944134e-16
cost at iteration 377 :  1.7210459237594076e-16
cost at iteration 378 :  1.7209635535175837e-16
cost at iteration 379 :  1.7137672687480183e-16
cost at iteration 380 :  1.7137910341159073e-16
cost at iteration 381 :  1.7066513222180986e-16
cost at iteration 382 :  1.706778428993043e-16
cost at iteration 383 :  1.699693141203524e-16
cost at iteration 384 :  1.699920889355304e-16
cost at iteration 385 :  1.6928879673878522e-16
cost at iteration 386 :  1.6932137467487976e-16
cost at iteration 387 :  1.6862312188564432e-16
cost at iteration 388 :  1.6866525050412874e-16
cost at iteration 389 :  1.6797184822807972e-16
cost at iteration 390 :  1.680232832820546e-16
cost at iteration 391 :  1.6733455055065764e-16
cost at iteration 392 :  1.6739505561858717e-16
cost at iteration 393 :  1.6671081905231273e-16
cost at iteration 394 :  1.6678016519078775e-16
cost at iteration 395 :  1.6610025867914903e-16
cost at iteration 396 :  1.6617822409369363e-16
cost at iteration 397 :  1.6550248849086405e-16
cost at iteration 398 :  1.6558885822372232e-16
cost at iteration 399 :  1.649171410591977e-16
cost at iteration 400 :  1.6501170669309652e-16
cost at iteration 401 :  1.6434386189606726e-16
cost at iteration 402 :  1.644464212731847e-16
cost at iteration 403 :  1.637823089102414e-16
cost at iteration 404 :  1.6389266586536326e-16
cost at iteration 405 :  1.6323215189041745e-16
cost at iteration 406 :  1.6335011599765204e-16
cost at iteration 407 :  1.6269307201367388e-16
cost at iteration 408 :  1.6281845834588197e-16
cost at iteration 409 :  1.621647613775169e-16
cost at iteration 410 :  1.6229739027774668e-16
cost at iteration 411 :  1.616469225544755e-16
cost at iteration 412 :  1.6178661941882892e-16
cost at iteration 413 :  1.611392681678837e-16
cost at iteration 414 :  1.6128586323910014e-16
cost at iteration 415 :  1.606415204876609e-16
cost at iteration 416 :  1.6079484865887437e-16
cost at iteration 417 :  1.6015341104511591e-16
cost at iteration 418 :  1.6031331167329163e-16
cost at iteration 419 :  1.596746802656835e-16
cost at iteration 420 :  1.598409969941008e-16
cost at iteration 421 :  1.5920507711866923e-16
cost at iteration 422 :  1.5937765770805937e-16
cost at iteration 423 :  1.5874435878306582e-16
cost at iteration 424 :  1.589230549508373e-16
cost at iteration 425 :  1.5829229032858286e-16
cost at iteration 426 :  1.5847695759587763e-16
cost at iteration 427 :  1.5784864441115632e-16
cost at iteration 428 :  1.5803914195717128e-16
cost at iteration 429 :  1.57413200982215e-16
cost at iteration 430 :  1.5760939150541144e-16
cost at iteration 431 :  1.5698574701080476e-16
cost at iteration 432 :  1.5718749659670622e-16
cost at iteration 433 :  1.565660762181502e-16
cost at iteration 434 :  1.5677325421340915e-16
cost at iteration 435 :  1.561539888238974e-16
cost at iteration 436 :  1.5636646771618433e-16
cost at iteration 437 :  1.5574929130339111e-16
cost at iteration 438 :  1.5596694660691191e-16
cost at iteration 439 :  1.5535179615563268e-16
cost at iteration 440 :  1.5557450630195892e-16
cost at iteration 441 :  1.5496132168116996e-16
cost at iteration 442 :  1.55188967915013e-16
cost at iteration 443 :  1.545776917695315e-16
cost at iteration 444 :  1.5481015804939887e-16
cost at iteration 445 :  1.5420073569578073e-16
cost at iteration 446 :  1.5443790859902958e-16
cost at iteration 447 :  1.5383028792560656e-16
cost at iteration 448 :  1.540720565579317e-16
cost at iteration 449 :  1.53466187928712e-16
cost at iteration 450 :  1.5371244383754085e-16
cost at iteration 451 :  1.5310827999991574e-16
cost at iteration 452 :  1.5335891709182028e-16
cost at iteration 453 :  1.527564130877784e-16
cost at iteration 454 :  1.5301132754944943e-16
cost at iteration 455 :  1.524104406302455e-16
cost at iteration 456 :  1.5266953085297048e-16
cost at iteration 457 :  1.5207022039702934e-16
cost at iteration 458 :  1.5233338690458242e-16
cost at iteration 459 :  1.5173561433851837e-16
cost at iteration 460 :  1.5200275971803373e-16
cost at iteration 461 :  1.5140648844066614e-16
cost at iteration 462 :  1.516775172766264e-16
cost at iteration 463 :  1.5108271258589096e-16
cost at iteration 464 :  1.5135753139687952e-16
cost at iteration 465 :  1.5076416041936353e-16
cost at iteration 466 :  1.5104267759757492e-16
cost at iteration 467 :  1.5045070922084224e-16
cost at iteration 468 :  1.5073283497409743e-16
cost at iteration 469 :  1.5014223978133123e-16
cost at iteration 470 :  1.5042788607760608e-16
cost at iteration 471 :  1.4983863628481296e-16
cost at iteration 472 :  1.501277167990092e-16
cost at iteration 473 :  1.4953978619449397e-16
cost at iteration 474 :  1.4983221625748982e-16
cost at iteration 475 :  1.4924558014345087e-16
cost at iteration 476 :  1.4954127669321576e-16
cost at iteration 477 :  1.4895591182962204e-16
cost at iteration 478 :  1.4925479336429774e-16
cost at iteration 479 :  1.4867067791469916e-16
cost at iteration 480 :  1.4897266444769876e-16
cost at iteration 481 :  1.4838977792698903e-16
cost at iteration 482 :  1.4869479094375846e-16
cost at iteration 483 :  1.481131141678875e-16
cost at iteration 484 :  1.4842107658454463e-16
cost at iteration 485 :  1.4784059162198223e-16
cost at iteration 486 :  1.4815142774546313e-16
cost at iteration 487 :  1.4757211787041374e-16
cost at iteration 488 :  1.4788575336029933e-16
cost at iteration 489 :  1.4730760300760668e-16
cost at iteration 490 :  1.4762396483936073e-16
cost at iteration 491 :  1.47046959560976e-16
cost at iteration 492 :  1.4736597599065448e-16
cost at iteration 493 :  1.467901024136654e-16
cost at iteration 494 :  1.4711170294397777e-16
cost at iteration 495 :  1.465369487301466e-16
cost at iteration 496 :  1.4686106407775486e-16
cost at iteration 497 :  1.4628741788446219e-16
cost at iteration 498 :  1.4661397994863944e-16
cost at iteration 499 :  1.4604143139116283e-16
cost at iteration 500 :  1.4637037322347936e-16
cost at iteration 501 :  1.4579891283868717e-16
cost at iteration 502 :  1.4613016861391273e-16
cost at iteration 503 :  1.455597878251272e-16
cost at iteration 504 :  1.458932928132044e-16
cost at iteration 505 :  1.45323983896376e-16
cost at iteration 506 :  1.456596744353838e-16
cost at iteration 507 :  1.4509143048637867e-16
cost at iteration 508 :  1.4542924395647864e-16
cost at iteration 509 :  1.4486205885955271e-16
cost at iteration 510 :  1.4520193365789582e-16
cost at iteration 511 :  1.4463580205525197e-16
cost at iteration 512 :  1.4497767757173801e-16
cost at iteration 513 :  1.4441259483410705e-16
cost at iteration 514 :  1.4475641142802497e-16
cost at iteration 515 :  1.4419237362632653e-16
cost at iteration 516 :  1.4453807260382347e-16
cost at iteration 517 :  1.4397507648170726e-16
cost at iteration 518 :  1.4432260007397844e-16
cost at iteration 519 :  1.4376064302145237e-16
cost at iteration 520 :  1.4410993436374158e-16
cost at iteration 521 :  1.4354901439157366e-16
cost at iteration 522 :  1.4390001750283616e-16
cost at iteration 523 :  1.4334013321792994e-16
cost at iteration 524 :  1.4369279298121054e-16
cost at iteration 525 :  1.431339435627649e-16
cost at iteration 526 :  1.4348820570620588e-16
cost at iteration 527 :  1.429303908827552e-16
cost at iteration 528 :  1.4328620196122045e-16
cost at iteration 529 :  1.4272942198842682e-16
cost at iteration 530 :  1.430867293657746e-16
cost at iteration 531 :  1.4253098500493843e-16
cost at iteration 532 :  1.4288973683682302e-16
cost at iteration 533 :  1.4233502933424897e-16
cost at iteration 534 :  1.4269517455142133e-16
cost at iteration 535 :  1.4214150561839265e-16
cost at iteration 536 :  1.4250299391064337e-16
cost at iteration 537 :  1.4195036570410074e-16
cost at iteration 538 :  1.4231314750464212e-16
cost at iteration 539 :  1.4176156260855069e-16
cost at iteration 540 :  1.4212558907882763e-16
cost at iteration 541 :  1.4157505048615555e-16
cost at iteration 542 :  1.4194027350125768e-16
cost at iteration 543 :  1.4139078459658481e-16
cost at iteration 544 :  1.4175715673092965e-16
cost at iteration 545 :  1.412087212736892e-16
cost at iteration 546 :  1.4157619578725473e-16
cost at iteration 547 :  1.410288178954645e-16
cost at iteration 548 :  1.4139734872042035e-16
cost at iteration 549 :  1.4085103285502554e-16
cost at iteration 550 :  1.412205745827154e-16
cost at iteration 551 :  1.4067532553248807e-16
cost at iteration 552 :  1.4104583340072967e-16
cost at iteration 553 :  1.4050165626764143e-16
cost at iteration 554 :  1.408730861485722e-16
cost at iteration 555 :  1.4032998633370535e-16
cost at iteration 556 :  1.4070229472178032e-16
cost at iteration 557 :  1.401602779116803e-16
cost at iteration 558 :  1.4053342191212164e-16
cost at iteration 559 :  1.399924940656708e-16
cost at iteration 560 :  1.4036643138316976e-16
cost at iteration 561 :  1.398265987188784e-16
cost at iteration 562 :  1.4020128764663507e-16
cost at iteration 563 :  1.3966255663039752e-16
cost at iteration 564 :  1.4003795603949572e-16
cost at iteration 565 :  1.3950033337271977e-16
cost at iteration 566 :  1.3987640270169748e-16
cost at iteration 567 :  1.3933989530991096e-16
cost at iteration 568 :  1.3971659455467772e-16
cost at iteration 569 :  1.391812095764909e-16
cost at iteration 570 :  1.3955849928045728e-16
cost at iteration 571 :  1.3902424405688303e-16
cost at iteration 572 :  1.3940208530141564e-16
cost at iteration 573 :  1.3886896736564823e-16
cost at iteration 574 :  1.3924732176066803e-16
cost at iteration 575 :  1.3871534882811324e-16
cost at iteration 576 :  1.3909417850302663e-16
cost at iteration 577 :  1.3856335846175704e-16
cost at iteration 578 :  1.3894262605652236e-16
cost at iteration 579 :  1.38412966958033e-16
cost at iteration 580 :  1.3879263561453963e-16
cost at iteration 581 :  1.382641456648277e-16
cost at iteration 582 :  1.386441790184188e-16
cost at iteration 583 :  1.3811686656939737e-16
cost at iteration 584 :  1.3849722874061082e-16
cost at iteration 585 :  1.3797110228179945e-16
cost at iteration 586 :  1.383517578683031e-16
cost at iteration 587 :  1.3782682601882337e-16
cost at iteration 588 :  1.3820774008754965e-16
cost at iteration 589 :  1.3768401158844114e-16
cost at iteration 590 :  1.380651496679117e-16
cost at iteration 591 :  1.3754263337463832e-16
cost at iteration 592 :  1.37923961447428e-16
cost at iteration 593 :  1.3740266632272429e-16
cost at iteration 594 :  1.3778415081813202e-16
cost at iteration 595 :  1.3726408592508702e-16
cost at iteration 596 :  1.3764569371194857e-16
cost at iteration 597 :  1.371268682073487e-16
cost at iteration 598 :  1.375085665870151e-16
cost at iteration 599 :  1.3699098971486164e-16