~libadjoint/dolfin-adjoint/web : revision 171

1

.. py:currentmodule:: dolfin_adjoint

2

3

==================================================

4

Mathematical Programs with Equilibrium Constraints

5

==================================================

6

7

.. sectionauthor:: Simon W. Funke <simon@simula.no>

8

9

This demo solves example X of :cite:`hintermueller2011`.

10

11

******************

12

Problem definition

13

******************

14

15

This problem is to minimise the dissipated power in the fluid

16

17

.. math::

18

\frac{1}{2} \int_{\Omega} \alpha(\rho) u \cdot u + \mu + \int_{\Omega} \nabla u : \nabla u - \int_{\Omega} f u

19

20

subject to the Stokes equations with velocity Dirichlet conditions

21

22

.. math::

23

\alpha(\rho) u - \mu \nabla^2 u + \nabla p &= f \qquad \mathrm{in} \ \Omega \\

24

\mathrm{div}(u) &= 0 \qquad \mathrm{on} \ \Omega \\

25

u &= b \qquad \mathrm{on} \ \delta \Omega \\

26

27

and to the control constraints on available fluid volume

28

29

.. math::

30

0 \le \rho(x) &\le 1 \qquad \forall x \in \Omega \\

31

\int_{\Omega} \rho &\le V

32

33

where :math:`u` is the velocity, :math:`p` is the pressure, :math:`\rho` is the control

34

(:math:`\rho(x) = 1` means fluid present, :math:`\rho(x) = 0` means no fluid present),

35

:math:`f` is a prescribed source term (here 0), :math:`V` is the volume bound on the control,

36

:math:`\alpha(\rho)` models the inverse permeability as a function

37

of the control

38

39

.. math::

40

\alpha(\rho) = \bar{\alpha} + (\underline{\alpha} - \bar{\alpha}) \rho \frac{1 + q}{\rho + q}

41

42

with :math:`\bar{\alpha}`, :math:`\underline{\alpha}` and :math:`q` prescribed

43

constants. The parameter :math:`q` penalises deviations from the values 0 or 1;

44

the higher q, the closer the solution will be to having the two discrete values 0 or 1.

45

46

The problem domain :math:`\Omega` is parameterised by the aspect ratio

47

:math:`\delta` (the domain is 1 unit high and :math:`\delta` units wide);

48

in this example, we will solve the harder problem of :math:`\delta = 1.5`.

49

The boundary conditions are specified in figure 10 of Borrvall and Petersson,

50

reproduced here.

51

52

.. image:: 2-stokes-topology-bcs.png

53

:scale: 80

54

:align: center

55

56

Physically, this problem corresponds to finding the fluid-solid distribution

57

:math:`\rho(x)` that minimises the dissipated power in the fluid.

58

59

As Borrvall and Petersson comment, it is necessary to solve this problem with

60

:math:`q=0.1` to ensure that the result approaches a discrete-valued solution,

61

but solving this problem directly with this value of :math:`q` leads to a local

62

minimum configuration of two straight pipes across the domain (like the top half

63

of figure 11). Therefore, we follow their suggestion to first solve the

64

optimisation problem with a smaller penalty parameter of :math:`q=0.01`; this

65

optimisation problem does not yield bang-bang solutions but is easier to solve,

66

and gives an initial guess from which the :math:`q=0.1` case converges to the

67

better minimum.

68

69

**************

70

Implementation

71

**************

72

73

First, the :py:mod:`dolfin` and :py:mod:`dolfin_adjoint` modules are imported:

74

75

.. code-block:: python

76

77

from dolfin import *

78

from dolfin_adjoint import *

79

80

Next we import the Python interface to IPOPT. If IPOPT is unavailable on your

81

system, we strongly :doc:`suggest you install it <../download/index>`; IPOPT is the leading open-source

82

optimisation algorithm.

83

84

.. code-block:: python

85

86

import pyipopt

87

88

Next we define some constants, and define the inverse permeability as a function

89

of :math:`\rho`.

90

91

.. code-block:: python

92

93

mu = Constant(1.0) # viscosity

94

alphaunderbar = 2.5 * mu / (100**2) # parameter for \alpha

95

alphabar = 2.5 * mu / (0.01**2) # parameter for \alpha

96

q = Constant(0.01) # q value that controls difficulty/

97

# discrete-valuedness of solution

98

99

def alpha(rho):

100

"""Inverse permeability as a function of rho, equation (40)"""

101

return alphabar + (alphaunderbar - alphabar) * rho * (1 + q) / (rho + q)

102

103

Next we define the mesh (a rectangle 1 high and :math:`\delta` wide)

104

and the function spaces to be used for the control :math:`\rho`,

105

the velocity :math:`u` and the pressure :math:`p`. Here we will use

106

the Taylor-Hood finite element to discretise the Stokes equations

107

:cite:`taylor1973`.

108

109

.. code-block:: python

110

111

n = 100

112

delta = 1.5 # The aspect ratio of the domain, 1 high and \delta wide

113

V = Constant(1.0/3) * delta # want the fluid to occupy 1/3 of the domain

114

115

mesh = RectangleMesh(0.0, 0.0, delta, 1.0, n, n)

116

A = FunctionSpace(mesh, "DG", 0) # control function space

117

U = VectorFunctionSpace(mesh, "CG", 2) # velocity function space

118

P = FunctionSpace(mesh, "CG", 1) # pressure function space

119

W = MixedFunctionSpace([U, P]) # mixed Taylor-Hood function space

120

121

Next we define the forward boundary condition.

122

123

.. code-block:: python

124

125

class InflowOutflow(Expression):

126

def eval(self, values, x):

127

values[1] = 0.0

128

values[0] = 0.0

129

l = 1.0/6.0

130

gbar = 1.0

131

132

if x[0] == 0.0 or x[0] == delta:

133

if (1.0/4 - l/2) < x[1] < (1.0/4 + l/2):

134

t = x[1] - 1.0/4

135

values[0] = gbar*(1 - (2*t/l)**2)

136

if (3.0/4 - l/2) < x[1] < (3.0/4 + l/2):

137

t = x[1] - 3.0/4

138

values[0] = gbar*(1 - (2*t/l)**2)

139

140

def value_shape(self):

141

return (2,)

142

143

Next we define a function that given a control :math:`\rho` solves the forward PDE

144

for velocity and pressure :math:`(u, p)`. (The advantage of formulating it in this manner

145

is that it makes it easy to conduct :doc:`Taylor remainder convergence tests

146

<../documentation/verification>`.)

147

148

.. code-block:: python

149

150

def forward(rho):

151

"""Solve the forward problem for a given fluid distribution rho(x)."""

152

w = Function(W)

153

(u, p) = split(w)

154

(v, q) = TestFunctions(W)

155

156

F = (alpha(rho) * inner(u, v) * dx + inner(grad(u), grad(v)) * dx +

157

inner(grad(p), v) * dx + inner(div(u), q) * dx)

158

bc = DirichletBC(W.sub(0), InflowOutflow(), "on_boundary")

159

solve(F == 0, w, bcs=bc)

160

161

return w

162

163

164

Now we define the ``__main__`` section. We define the initial guess for the

165

control and use it to solve the forward PDE. In order to ensure feasibility of

166

the initial control guess, we interpolate the volume bound; this ensures that

167

the integral constraint and the bound constraint are satisfied.

168

169

.. code-block:: python

170

171

if __name__ == "__main__":

172

rho = interpolate(Constant(float(V)/delta), A, name="Control")

173

w = forward(rho)

174

(u, p) = split(w)

175

176

With the forward problem solved once, :py:mod:`dolfin_adjoint` has built a

177

*tape* of the forward model; it will use this tape to drive the optimisation, by

178

repeatedly solving the forward model and the adjoint model for varying control

179

inputs.

180

181

As in the :doc:`Poisson topology example <1-poisson-topology>`, we will use an

182

evaluation callback to dump the control iterates to disk for visualisation. As

183

this optimisation problem (:math:`q=0.01`) is solved only to generate an initial

184

guess for the main task (:math:`q=0.1`), we shall save these iterates in

185

``output/control_iterations_guess.pvd``.

186

187

.. code-block:: python

188

189

controls = File("output/control_iterations_guess.pvd")

190

rho_viz = Function(A, name="ControlVisualisation")

191

def eval_cb(j, rho):

192

rho_viz.assign(rho)

193

controls << rho_viz

194

195

Now we define the functional and :doc:`reduced functional <../maths/2-problem>`:

196

197

.. code-block:: python

198

199

J = Functional(0.5 * inner(alpha(rho) * u, u) * dx + mu * inner(grad(u), grad(u)) * dx)

200

m = SteadyParameter(rho)

201

Jhat = ReducedFunctional(J, m, eval_cb=eval_cb)

202

rfn = ReducedFunctionalNumPy(Jhat)

203

204

The control constraints are the same as the :doc:`Poisson topology example

205

<1-poisson-topology>`, and so won't be discussed again here.

206

207

.. code-block:: python

208

209

# Bound constraints

210

lb = 0.0

211

ub = 1.0

212

213

class VolumeConstraint(InequalityConstraint):

214

"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""

215

def __init__(self, V):

216

self.V = float(V)

217

218

# The derivative of the constraint g(x) is constant

219

# (it is the negative of the diagonal of the lumped mass matrix for the

220

# control function space), so let's assemble it here once.

221

# This is also useful in rapidly calculating the integral each time

222

# without re-assembling.

223

self.smass = assemble(TestFunction(A) * Constant(1) * dx)

224

self.tmpvec = Function(A)

225

226

def function(self, m):

227

self.tmpvec.vector()[:] = m

228

229

# Compute the integral of the control over the domain

230

integral = self.smass.inner(self.tmpvec.vector())

231

print "Current control: ", integral

232

return [self.V - integral]

233

234

def jacobian(self, m):

235

return [-self.smass]

236

237

def length(self):

238

"""Return the number of components in the constraint vector (here, one)."""

239

return 1

240

241

Now that all the ingredients are in place, we can perform the initial

242

optimisation. We set the maximum number of iterations for this initial

243

optimisation problem to 30; there's no need to solve this to completion,

244

as its only purpose is to generate an initial guess.

245

246

.. code-block:: python

247

248

# Solve the optimisation problem with q = 0.01

249

nlp = rfn.pyipopt_problem(bounds=(lb, ub), constraints=VolumeConstraint(V))

250

nlp.int_option('max_iter', 30)

251

rho_opt = nlp.solve()

252

File("output/control_solution_guess.xml.gz") << rho_opt

253

254

With the optimised value for :math:`q=0.01` in hand, we *reset* the

255

dolfin-adjoint state, clearing its tape, and configure the new problem

256

we want to solve. We need to update the values of :math:`q` and

257

:math:`\rho`:

258

259

.. code-block:: python

260

261

q.assign(0.1)

262

rho.assign(rho_opt)

263

adj_reset()

264

265

Since we have cleared the tape, we need to execute the forward model once again

266

to redefine the problem. (It is also possible to modify the tape, but this way

267

is easier to understand.) We will also redefine the functionals and parameters;

268

this time, the evaluation callback will save the optimisation iterations to

269

``output/control_iterations_final.pvd``.

270

271

.. code-block:: python

272

273

w = forward(rho)

274

(u, p) = split(w)

275

276

# Define the reduced functionals

277

controls = File("output/control_iterations_final.pvd")

278

rho_viz = Function(A, name="ControlVisualisation")

279

def eval_cb(j, rho):

280

rho_viz.assign(rho)

281

controls << rho_viz

282

283

J = Functional(0.5 * inner(alpha(rho) * u, u) * dx + mu * inner(grad(u), grad(u)) * dx)

284

m = SteadyParameter(rho)

285

Jhat = ReducedFunctional(J, m, eval_cb=eval_cb)

286

rfn = ReducedFunctionalNumPy(Jhat)

287

288

We can now solve the optimisation problem with :math:`q=0.1`, starting from the

289

solution of :math:`q=0.01`:

290

291

.. code-block:: python

292

293

nlp = rfn.pyipopt_problem(bounds=(lb, ub), constraints=VolumeConstraint(V))

294

nlp.int_option('max_iter', 100)

295

rho_opt = nlp.solve()

296

File("output/control_solution_final.xml.gz") << rho_opt

297

298

The example code can be found in ``examples/stokes-topology/`` in the ``dolfin-adjoint`` source tree,

299

and executed as follows:

300

301

.. code-block:: bash

302

303

$ mpiexec -n 4 python stokes-topology.py

304

...

305

Number of Iterations....: 100

306

307

(scaled) (unscaled)

308

Objective...............: 4.5944633030224409e+01 4.5944633030224409e+01

309

Dual infeasibility......: 1.8048641504211900e-03 1.8048641504211900e-03

310

Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00

311

Complementarity.........: 9.6698653740681504e-05 9.6698653740681504e-05

312

Overall NLP error.......: 1.8048641504211900e-03 1.8048641504211900e-03

313

314

315

Number of objective function evaluations = 105

316

Number of objective gradient evaluations = 101

317

Number of equality constraint evaluations = 0

318

Number of inequality constraint evaluations = 105

319

Number of equality constraint Jacobian evaluations = 0

320

Number of inequality constraint Jacobian evaluations = 101

321

Number of Lagrangian Hessian evaluations = 0

322

Total CPU secs in IPOPT (w/o function evaluations) = 11.585

323

Total CPU secs in NLP function evaluations = 556.795

324

325

EXIT: Maximum Number of Iterations Exceeded.

326

327

The optimisation iterations can be visualised by opening

328

``output/control_iterations_final.pvd`` in

329

paraview. The resulting solution appears very similar to the solution proposed

330

in :cite:`borrvall2003`.

331

332

.. image:: 2-stokes-topology.png

333

:scale: 25

334

:align: center

335

336

.. rubric:: References

337

338

.. bibliography:: /documentation/3-mpec.bib

339

:cited:

340

:labelprefix: 2E