Abstract
We present a convection–diffusion inverse problem that aims to identify an unknown number of sources and their locations. We model the sources using a binary function, and we show that the inverse problem can be formulated as a large-scale mixed-integer nonlinear optimization problem. We show empirically that current state-of-the-art mixed-integer solvers cannot solve this problem and that applying simple rounding heuristics to solutions of the relaxed problem can fail to identify the correct number and location of the sources. We develop two new rounding heuristics that exploit the value and a physical interpretation of the continuous relaxation solution, and we apply a steepest-descent improvement heuristic to obtain satisfactory solutions to both two- and three-dimensional inverse problems. We also provide the code used in our numerical experiments in open-source format.
Similar content being viewed by others
References
Abhishek K, Leyffer S, Linderoth JT (2010) FilMINT: an outer-approximation-based solver for nonlinear mixed integer programs. INFORMS J Comput 22:555–567. https://doi.org/10.1287/ijoc.1090.0373
Achterberg T (2005) SCIP—a framework to integrate constraint and mixed integer programming. Technical Report ZIB-Report 04-19, Konrad-Zuse-Zentrum für Informationstechnik Berlin, Takustr. 7, Berlin
Achterberg T (2009) Scip: solving constraint integer programs. Math Program Comput 1(1):1–41
Akcelik V, Biros G, Draganescu A, Ghattas O, Hill J, van Bloemen Waanders B (2005) Dynamic data-driven inversion for terascale simulations: Real-time identification of airborne contaminants. In: Proceedings of SC2005, Seattle, WA
Akçelik V, Biros G, Ghattas O, Hill J, Keyes D, van Bloemen Waanders B (2006) Parallel algorithms for PDE-constrained optimization. In: Parallel processing for scientific computing. SIAM, pp 291–322
Akrotirianakis I, Maros I, Rustem B (2001) An outer approximation based branch-and-cut algorithm for convex 0–1 MINLP problems. Optim Methods Softw 16:21–47
Ascher UM, Haber E (2001) Grid refinement and scaling for distributed parameter estimation problems. Inverse Prob 17:571–590
Balas E (1975) Facets of the knapsack polytope. Math Program 8:146–164
Bangerth W, Klie H, Matossian V, Parashar M, Wheeler MF (2005) An autonomic reservoir framework for the stochastic optimization of well placement. Cluster Comput 8(4):255–269
Bangerth W, Klie H, Wheeler M, Stoffa P, Sen M (2006) On optimization algorithms for the reservoir oil well placement problem. Comput Geosci 10(3):303–319. https://doi.org/10.1007/s10596-006-9025-7
Bartlett R, Heinkenschloss M, Ridzal D, van Bloemen Waanders B (2005) Domain decomposition methods for advection dominated linear-quadratic elliptic optimal control problems. Comput Methods Appl Mech Eng 195(44–47):6428–6447
Bellout MC, Ciaurri DE, Durlofsky LJ, Foss B, Kleppe J (2012) Joint optimization of oil well placement and controls. Comput Geosci 16(4):1061–1079
Belotti P, Kirches C, Leyffer S, Linderoth J, Luedtke J, Mahajan A (2013) Mixed-integer nonlinear optimization. Acta Numerica 22:1–131. https://doi.org/10.1017/S0962492913000032
Belotti P, Kirches C, Leyffer S, Linderoth J, Luedtke J, Mahajan A (2013) Mixed integer nonlinear programming. Acta Numerica 22:1–131
Bendsøe M, Sigmund O (2004) Topological optimization theory. Springer, Berlin
Bezanson J, Karpinski S, Shah VB, Edelman A (2012) Julia: a fast dynamic language for technical computing. arXiv preprint arXiv:1209.5145
Biegler L, Ghattas O, Heinkenschloss M, van Bloemen Waanders B (eds) (2001) Large-scale PDE-constrained optimization, vol 30. Lecture notes in computational science and engineering. Springer, Berlin
Biegler L, Ghattas O, Heinkenschloss M, Keyes D, van Bloemen Waanders B (eds) (2007) Real-time PDE-constrained optimization. SIAM, Philadelphia
Biegler LT, Ghattas O, Heinkenschloss M, van Bloemen Waanders B (2003) Large-scale PDE-constrained optimization: an introduction. Large-scale PDE-constrained optimization. Springer, Berlin, pp 3–13
Biros G, Ghattas O (2005) Parallel Lagrange-Newton-Krylov-Schur Methods for PDE-constrained optimization. Part I: The Krylov-Schur Solver. SIAM J Sci Comput 27(2):687–713
Biros G, Ghattas O (2005) Parallel Lagrange-Newton-Krylov-Schur Methods for PDE-Constrained Optimization, Part II: The Lagrange-Newton Solver, and its Application to Optimal Control of Steady Viscous Flows. SIAM J Sci Comput 27(2):714–739
Bonami P, Biegler L, Conn A, Cornuéjols G, Grossmann I, Laird C, Lee J, Lodi A, Margot F, Sawaya N, Wächter A (2008) An algorithmic framework for convex mixed integer nonlinear programs. Discrete Optim 5(2):186–204
Bonami P, Cornuéjols G, Lodi A, Margot F (2009) A feasibility pump for mixed integer nonlinear programs. Math Program 119:331–352
Bonami P, Lee J (2007) Bonmin user’s manual. Numer Math 4:1–32
Borzi A (2007) High-order discretization and multigrid solution of elliptic nonlinear constrained optimal control problems. J Comp Appl Math 200:67–85
Borzì A, Schulz V (2009) Multigrid methods for PDE optimization. SIAM Rev 51(2):361–395. https://doi.org/10.1137/060671590
Burer S, Letchford A (2012) Non-convex mixed-integer nonlinear programming: a survey. Surv Oper Res Manag Sci 17:97–106
Bürger A, Zeile C, Hahn M, Altmann-Dieses A, Sager S, Diehl M (2020) pycombina: an open-source tool for solving combinatorial approximation problems arising in mixed-integer optimal control. In: IFAC World Congress 2020. Accepted
Bussieck MR, Pruessner A (2003) Mixed-integer nonlinear programming. SIAG/OPT Views-and-News 14(1):19–22
Çezik MT, Iyengar G (2005) Cuts for mixed 0–1 conic programming. Math Program 104:179–202
Chan TF, Shen J (2005) Image Processing and Analysis. Society for Industrial and Applied Mathematics. http://bookstore.siam.org/ot94/
Committee T (2010) Advanced fueld pellet materials and fuel rod design for water cooled reactors. Technical report, International Atomic Energy Agency
Conn AR, Gould NI, Toint PL (2000) Trust region methods, vol 1. SIAM, Philadelphia
Costa MFP, Rocha AMA, Francisco RB, Fernandes EM (2016) Firefly penalty-based algorithm for bound constrained mixed-integer nonlinear programming. Optimization 65(5):1085–1104
Csurka G, Larlus D, Perronnin F, Meylan F (2013) What is a good evaluation measure for semantic segmentation? In: BMVC, vol 27. Citeseer
Dakin RJ (1965) A tree search algorithm for mixed programming problems. Comput J 8:250–255
Danna E, Rothberg E, LePape C (2005) Exploring relaxation induced neighborhoods to improve MIP solutions. Math Program 102:71–90
De Wolf D, Smeers Y (2000) The gas transmission problem solved by an extension of the simplex algorithm. Manage Sci 46:1454–1465
Donovan G, Rideout D (2003) An integer programming model to optimize resource allocation for wildfire containment. Forest Sci 61(2):331–335
Drewes S (2009) Mixed integer second order cone programming. Ph.D. thesis, Technische Universität Darmstadt
Drewes S, Ulbrich S (2012) Subgradient based outer approximation for mixed integer second order cone programming. In: Mixed integer nonlinear programming, The IMA volumes in mathematics and its applications, vol 154. Springer, New York, pp 41–59. ISBN 978-1-4614-1926-6
Dunning I, Huchette J, Lubin M (2017) Jump: a modeling language for mathematical optimization. SIAM Rev 59(2):295–320
Duran MA, Grossmann I (1986) An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math Program 36:307–339
Ehrhardt K, Steinbach MC (2005) Nonlinear optimization in gas networks. Springer, Berlin
Engl HW, Hanke M, Neubauer A (1996) Regularization of inverse problems, mathematics and its applications, vol 375. Kluwer Academic Publishers Group, Dordrecht. https://doi.org/10.1007/978-94-009-1740-8
Fipki S, Celi A (2008) The use of multilateral well designs for improved recovery in heavy oil reservoirs. In: IADV/SPE Conference and Exhibition. SPE, Orlanda, Florida
Fischetti M, Glover F, Lodi A (2005) The feasibility pump. Math Program 104:91–104
Fischetti M, Lodi A (2002) Local branching. Math Program 98:23–47
Floudas CA (2000) Deterministic global optimization: theory, algorithms and applications. Kluwer Academic Publishers, Berlin
Fourer R, Gay DM, Kernighan BW (1993) AMPL: a modeling language for mathematical programming. The Scientific Press, Beijing
Frangioni A, Gentile C (2006) Perspective cuts for a class of convex 0–1 mixed integer programs. Math Program 106:225–236
Fügenschuh A, Geißler B, Martin A, Morsi A (2009) The transport PDE and mixed-integer linear programming. In: Dagstuhl Seminar Proceedings. Schloss Dagstuhl-Leibniz-Zentrum für Informatik
Garmatter D, Porcelli M, Rinaldi F, Stoll M (2019) Improved penalty algorithm for mixed integer PDE constrained optimization (MIPDECO) problems. arXiv preprint arXiv:1907.06462
Geoffrion AM (1972) Generalized Benders decomposition. J Optim Theory Appl 10(4):237–260
Gerdts M, Sager S (2012) Mixed-integer DAE optimal control problems: necessary conditions and bounds. In: Biegler L, Campbell S, Mehrmann V (eds) Control and optimization with differential-algebraic constraints. SIAM, Berlin, pp 189–212
Golub GH, Heath M, Wahba G (1979) Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21(2):215–223. https://doi.org/10.1080/00401706.1979.10489751
Grossmann IE (2002) Review of nonlinear mixed-integer and disjunctive programming techniques. Optim Eng 3:227–252
Grossmann IE, Kravanja Z (1997) Mixed-integer nonlinear programming: a survey of algorithms and applications. In: Biegler ACLT, Coleman TF, Santosa F (eds) Large-scale optimization with applications, Part II: optimal design and control. Springer, New York
Gugat M, Leugering G, Martin A, Schmidt M, Sirvent M, Wintergerst D (2018) MIP-based instantaneous control of mixed-integer PDE-constrained gas transport problems. Comput Optim Appl 70(1):267–294
Günlük O, Linderoth J (2008) Perspective relaxation of mixed integer nonlinear programs with indicator variables. In: Lodi A, Panconesi A, Rinaldi G (eds.) IPCO 2008: The thirteenth conference on integer programming and combinatorial optimization, vol 5035, pp. 1–16
Gunzburger MD (2003) Perspectives in flow control and optimization, advances in design and control, vol 5. Society for Industrial and Applied Mathematics (SIAM), Philadelphia
Guo Jy Lu, Wx Yang Qc, Ts Miao (2019) The application of 0–1 mixed integer nonlinear programming optimization model based on a surrogate model to identify the groundwater pollution source. J Contam Hydrol 220:18–25
Gupta OK, Ravindran A (1985) Branch and bound experiments in convex nonlinear integer programming. Manage Sci 31:1533–1546
Haber E, Ascher UM (2001) Preconditioned all-at-once methods for large, sparse parameter estimation problems. Inverse Prob 17:1847–1864
Haber E, Oldenburg D (2000) A GCV based method for nonlinear ill-posed problems. Comput Geosci 4:41–63. https://doi.org/10.1023/A:1011599530422
Hahn M, Leyffer S, Zavala VM (2017) Mixed-Integer PDE-Constrained Optimal Control of Gas Networks. Technical Report Preprint ANL/MCS-P7095-0817, Mathematics and Computer Science Division, Argonne National Laboratory
Hahn M, Sager S, Leyffer S (2020) Binary optimal control by trust-region steepest descent. http://www.optimization-online.org/DB_FILE/2020/01/7589.pdf. Submitted for publication
Hansen PC (1998) Rank-deficient and discrete ill-posed problems. SIAM monographs on mathematical modeling and computation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia. https://doi.org/10.1137/1.9780898719697
Hante FM (2017) Relaxation methods for hyperbolic PDE mixed-integer optimal control problems. Opt Control Appl Methods 38(6):1103–1110
Hante FM, Sager S (2013) Relaxation methods for mixed-integer optimal control of partial differential equations. Comput Optim Appl 55(1):197–225
Hazra SB, Schulz V (2006) Simultaneous pseudo-timestepping for aerodynamic shape optimization problems with state constraints. SIAM J Sci Comput 28:1078–1099
Heinkenschloss M, Ridzal D (2008) Lecture notes in computational science and engineering, chapter. Integration of sequential quadratic programming and domain decomposition methods for nonlinear optimal control problems. Springer, Berlin
Herrmann M, Herzog R, Schmidt S, Vidal-Núñez J, Wachsmuth G (2018) Discrete total variation with finite elements and applications to imaging. arXiv.org
Hintermuller M, Vicente LN (2005) Space mapping for optimal control of partial differential equations. SIAM J Opt 15:1002–1025
Hinze M, Pinnau R, Ulbrich M, Ulbrich S (2009) Optimization with PDE constraints. Springer, Berlin
Horowitz E, Sahni S (1974) Computing partitions with applications to the knapsack problem. J ACM 21:277–292
Jeroslow RG (1973) There cannot be any algorithm for integer programming with quadratic constraints. Oper Res 21(1):221–224
Jung M (2013) Relaxations and approximations for mixed-integer optimal control. Ph.D. thesis, University Heidelberg. http://www.ub.uni-heidelberg.de/archiv/16036
Jung M, Reinelt G, Sager S (2015) The Lagrangian relaxation for the combinatorial integral approximation problem. Optim Methods Softw 30(1):54–80
Kannan R, Monma C (1978) On the computational complexity of integer programming problems. In: Henn R, Korte B, Oettli W (eds.) Optimization and Operations Research. Lecture notes in economics and mathematical systems, vol 157. Springer, Berlin, pp 161–172
Laird CD, Biegler LT, van Bloemen Waanders B, Bartlett RA (2005) Time dependent contaminant source determination for municipal water networks using large scale optimization. ASCE J Water Res Mgt Plan pp. 125–134
Lee J, Leyffer S (eds) (2011) Mixed integer nonlinear programming, IMA volume in mathematics and its applications. Springer, New York
Legg M, Davidson R, Nozick L (2013) Optimization-based regional hurricane mitigation planning. J Infrastruct Syst 19:1–11
Leyffer S, Munson T, Wild S, van Bloemen Waanders B, Ridzal D (2013) Mixed-integer PDE-constrained optimization. Position Paper #15 submitted in response to the ExaMath13 Call for Position Papers. https://collab.mcs.anl.gov/download/attachments/7569466/examath13_submission_15.pdf
Lucidi S, Rinaldi F (2013) An exact penalty global optimization approach for mixed-integer programming problems. Optim Lett 7(2):297–307
Mahajan A, Leyffer S, Linderoth J, Luedtke J, Munson T (2011) MINOTAUR: a toolkit for solving mixed-integer nonlinear optimization. wiki-page. http://wiki.mcs.anl.gov/minotaur
Mahajan A, Leyffer S, Linderoth J, Luedtke J, Munson T (2017) Minotaur: a mixed-integer nonlinear optimization toolkit. Technical reprt, ANL/MCS-P8010-0817, Argonne National Laboratory
Manns P, Kirches C (2018) Multi-dimensional sum-up rounding for elliptic control systems. DFG SPP 1962 Preprint . https://spp1962.wias-berlin.de/preprints/080.pdf. (submitted to SIAM Journal on Numerical Analysis)
Manns P, Kirches C (2019) Improved regularity assumptions for partial outer convexification of mixed-integer PDE-constrained optimization problems. ESAIM: control, optimisation and calculus of variations. http://www.optimization-online.org/DB_HTML/2018/04/6585.html. (accepted)
Manns P, Kirches C (2019) Multi-dimensional sum-up rounding using Hilbert curve iterates. In: Proceedings in applied mathematics and mechanics. (accepted)
Martello S, Pisinger D, Toth P (1999) Dynamic programming and strong bounds for the 0–1 knapsack problem. Manage Sci 45(3):414–424
Martello S, Toth P (1988) A new algorithm for the 0–1 knapsack problem. Manage Sci 34:633–644
Martello S, Toth P (1990) Knapsack problems: algorithms and computer implementations. Wiley, Chichester
Martin A, Möller M, Moritz S (2006) Mixed integer models for the stationary case of gas network optimization. Math Program 105:563–582
Nannicini G, Belotti P, Liberti L (2008) A local branching heuristic for MINLPs. arXiv:0812.2188v1 [math.CO]
Nocedal J, Wright S (2000) Numerical optimization. Springer, Berlin
Ozdogan U (2004) Optimization of well placement under time-dependent uncertainty. Master’s thesis, Stanford University
Pisinger D (1995) An expanding-core algorithm for the exact 0–1 knapsack problem. Eur J Oper Res 87:175–187
Pisinger D, Toth P (1998) Knapsack problems. In: Du DZ, Pardalos P (eds) Handbook of combinatorial optimization. Kluwer, Berlin, pp 1–89
Quesada I, Grossmann IE (1992) An LP/NLP based branch-and-bound algorithm for convex MINLP optimization problems. Comput Chem Eng 16:937–947
Reinke CM, la Mata Luque TMD, Su MF, Sinclair MB, El-Kady I (2011) Group-theory approach to tailored electromagnetic properties of metamaterials: an inverse-problem solution. Phys Rev E 83(6):06660–1–18 3
Rudin LI, Osher S, Fatemi E (1992) Nonlinear total variation based noise removal algorithms. Physica D 60(1–4):259–268. https://doi.org/10.1016/0167-2789(92)90242-F
Ruthotto L, Treister E, Haber E (2017) jinv-a flexible julia package for pde parameter estimation. SIAM J Sci Comput 39(5):S702–S722
Sager S (2006) Numerical methods for mixed–integer optimal control problems. Ph.D. thesis, Universität Heidelberg. https://mathopt.de/PUBLICATIONS/Sager2005.pdf
Sager S (2009) Reformulations and algorithms for the optimization of switching decisions in nonlinear optimal control. J Process Control 19(8):1238–1247 https://mathopt.de/PUBLICATIONS/Sager2009b.pdf
Sager S, Bock H, Diehl M (2012) The integer approximation error in mixed-integer optimal control. Math Program A 133(1–2):1–23 https://mathopt.de/PUBLICATIONS/Sager2012a.pdf
Sager S, Jung M, Kirches C (2011) Combinatorial integral approximation. Math Methods Oper Res 73(3):363–380. https://doi.org/10.1007/s00186-011-0355-4
Scherzer O, Grasmair M, Grossauer H, Haltmeier M, Lenzen F (2013) Variational methods in imaging. Springer, Berlin
Sharma S (2013) Mixed-integer nonlinear programming heuristics applied to a shale gas production optimization problem. Master’s thesis, Norwegian University of Science and Technology. http://www.diva-portal.org/smash/get/diva2:646797/FULLTEXT01.pdf
Sigmund O, Maute K (2013) Topological optimization approaches. Struct Multidiscip Opt 48:1031–1055
Sigmund O, Maute K (2013) Topology optimization approaches: a comparative review. Struct Multidiscip Optim 48(6):1031–1055
Simon R (2008) Multigrid solver for saddle point problems in PDE-constrained optimization. Ph.D. thesis, Johannes Kepler Universitat Linz
Steinbach MC (2007) On PDE solution in transient optimization of gas networks. J Comput Appl Math 203(2):345–361
Still C, Westerlund T (2006) Solving convex MINLP optimization problems using a sequential cutting plane algorithm. Comput Optim Appl 34(1):63–83
Stubbs R, Mehrotra S (1999) A branch-and-cut method for 0–1 mixed convex programming. Math Program 86:515–532
Tawarmalani M, Sahinidis NV (2002) Convexification and global optimization in continuous and mixed-integer nonlinear programming: theory, algorithms, software, and applications. Kluwer Academic Publishers, Boston
Vogel CR (1999) Sparse matrix computations arising in distributed parameter identification. SIAM J Matrix Anal Appl 20:1027–1037
Vogel CR (2002) Computational methods for inverse problems. Soc Ind Appl Math. doi 10(1137/1):9780898717570
You F, Leyffer S (2010) Oil spill response planning with MINLP. SIAG/OPT Views-and-News 21(2):1–8
You F, Leyffer S (2011) Mixed-integer dynamic optimization for oil-spill response planning with integration of a dynamic oil weathering model. AIChe J. https://doi.org/10.1002/aic.12536
Yu J, Anitescu M (2019) Multidimensional sum-up rounding for integer programming in optimal experimental design. Mathe. https://doi.org/10.1007/s10107-019-01421-z
Zavala VM (2014) Stochastic optimal control model for natural gas networks. Comput Chem Eng 64:103–113
Zhang P, Romero D, Beck J, Amon C (2013) Integration of AI and OR techniques in constraint programming for combinatorial optimization problems, chapter solving wind farm layout optimization with mixed integer programming and constraint programming. Springer, Berlin
Acknowledgements
This work was initiated as a part of the SAMSI Program on Optimization 20162017. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. Part of this material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under Contract DE-AC02-06CH11357. This work was also supported by the U.S. Department of Energy through grant DE-FG02-05ER25694. Ruthotto’s work was partially supported by the US National Science Foundation award DMS 1522599. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. SAND2019-10626 J. Hahn’s work was supported by the DFG (German Research Foundation) under SPP 1962 and by the German Federal Ministry of Education and Research (BMBF) under P2Chem grant 05M18NMA. We are grateful to Prof. Martin Siebenborn who helped us with the FEM derivation and deriving the weak form of the PDE. The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-access-plan
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1: Finite-element discretization of the source inversion problem
Here, we briefly review how we discretize the variational problem (1) on a computational mesh. We employ the first-discretize-then-optimize approach, which is a common strategy in PDE-constrained optimization; see; for example, (Gunzburger 2003, Sec. 2.9).
We begin by discretizing the PDE-constraint in (1) following the usual procedure of finite-element methods. To obtain a weak form of the constraint, we multiply both sides of the PDE constraint in (1) with a test function \(\phi \in H^1(\varOmega ,\mathbb R)\), the Sobolev space of all functions whose first derivative is square integrable over \(\varOmega\), with \(\phi (x) = 0\) for \(x\in \varGamma _D\), \(\frac{\partial \phi (x)}{\partial n} = 0\) for \(x\in \varGamma _N\), and then apply Green’s identity and integration by parts:
where the last identity is obtained by using the Neumann boundary conditions on \(\varGamma _N\) and the fact that \(\phi (x)=0\) on \(\varGamma _D\). The weak problem then is to find a \(u \in H^1(\varOmega ,\mathbb R)\) that satisfies (19) for all test functions \(\phi \in H^1(\varOmega ,\mathbb R)\).
Next, we discretize the state u and the control w in the weak form of the PDE constraint (19) using compactly supported Ansatz functions on a computational mesh. For ease of presentation we assume that our computational domain \(\varOmega = [0,1]^d\) is divided into \(N^d\) quadrilateral finite elements of edge length \(L=1/N\). We assume a lexicographical ordering of the elements \(\varOmega _1, \ldots , \varOmega _{N^d}\) in the mesh, which in our case consists of pixels and voxels for \(d=2,3\), respectively. We note that extensions to more general domains and anisotropic or unstructured meshes are straightforward. For example, the implementation used in our numerical experiments uses rectangular meshes whose elements have different edge lengths along the coordinate direction. We use the standard bi-/trilinear Ansatz functions \(\phi _1, \phi _2, \ldots ,\phi _{(N+1)^d}\) for \(d=2\) and \(d=3\), respectively, as test functions and to discretize the state variable u:
Because the Ansatz functions form a Lagrange basis (i.e., for the jth node of the mesh we have \(\phi _i(x_j) = \delta _{ij}\)), the elements in \(\mathbf {u}\) correspond to the value of u at the nodes of the mesh. We represent the control variable w as a linear combination of piecewise constant Ansatz functions \(\psi _1, \psi _2, \ldots , \psi _{N^d}\):
Similar to above, the Ansatz functions form a Lagrange basis (i.e., \(\psi _i(x)=1\) if \(x\in \varOmega _i\) and \(\psi (x)=0\) else); hence the entries in \(\mathbf {w}\) correspond to the values of w in the pixels/voxels of our mesh.
The finite-dimensional approximations of u and w in (20) and (21) allow us to write the weak form of the PDE-constraint (19) in terms of the coefficient vectors \(\mathbf {u}\) and \(\mathbf {w}\) as
where \(\mathbf {S}\in \mathbb R^{(N+1)^d \times (N+1)^d}\) is the (nonsingular) stiffness and \(\mathbf {M}\in \mathbb R^{(N+1)^d \times N^d}\) is the (full-rank) mass matrix, respectively. We compute the entries of both matrices, \(\mathbf {M}\) and \(\mathbf {S}\), approximately by applying a 5th-order Gaussian quadrature rule to (19):
Because of the compact support of the basis functions, the integrands vanish on all but a few elements, which leads to both matrices being sparse.
We now derive our discretization of the total variation regularizer (3). In this work, we use a first-order finite-difference discretization of the regularizer. This choice is motivated by its simplicity and computational efficiency. We refer to the recent work in Herrmann et al. (2018) for a more extensive discussion and finite-element discretizations of the total variation. Because we assume a quadrilateral mesh, the discrete gradient operators for \(d=2\) and \(d=3\) are
Here, \(\mathbf {I}_N\in \mathbb R^{N\times N}\) is the identity matrix, \(\otimes\) denotes the Kronecker product, and the one-dimensional difference operator is
We note that this discretization involves a grid change (from the N cell-centers to the \(N+1\) nodes). To accommodate for the grid change, the relaxed form includes the average matrix \(\mathbf {A}\)
for \(d=2\) and
for \(d=3\), respectively, where we use the one-dimensional average matrix
This leads to the discrete regularization function
where \(\mathbf {e}\) is a vector of all ones, the factor \(L^d\) represents the volume of the cells, and \(\kappa >0\) is a conditioning parameter (in our experiments, we use \(\kappa =10^{-3}\)). Note that the square and the square root are applied component-wise.
Appendix 2: Finite-difference discretization of the source inversion problem
We discretize both the source, w, and the state, u, in the cell-centered points of our \(N_x \times N_y\) computational mesh:
where \(L_x=2/N_x\) and \(L_y=1/N_y\) are the discretization steps in the x and y direction, respectively, and \(k=1,\ldots ,N_x, l=1,\ldots ,N_y\), \(i=0,\ldots ,N_x+1, j=0,\ldots ,N_y\). Here, the variables \(\mathbf {U}_{i,0}, \mathbf {U}_{N_x+1,j}, \mathbf {U}_{i,N_y+1}\) approximate the PDE solution at ghost points placed along \(\varGamma _N\), and the variables \(\mathbf {U}_{0,j}\) are the ghost points near \(\varGamma _D\).
Let us further define \(\mathbf {V}\in \mathbb R^{m}\) to obtain the bilinear interpolation from the cell-centered points of the mesh closest to the receiver locations \(r^1, r^2, \ldots , r^m\). Mathematically, for each receiver location \(r^k=(r^k_x, r^k _y), k=1,\ldots , m\), we define variable \(V_k\) as
where, \(i = \lfloor \dfrac{r^k_x}{L_x} + 0.5 \rfloor\) and \(j = \lfloor \dfrac{r^k_y}{L_y} + 0.5 \rfloor\), leading to the following mixed-integer quadratic program:
Here, the fifth row explicitly encodes the Neumann and Dirichlet boundary conditions. As before, we can again eliminate the state variables \(\mathbf {U}\) using the discretized PDE and boundary conditions and the state variables \(\mathbf {V}\) using (28), resulting in a problem that has similar structure to (6). As before, the objective function is second-order cone representable.
Appendix 3: Selection of regularization parameter for the relaxed problem
To find an effective regularization parameter, we consider the continuous relaxation (7) and follow the L-curve procedure; see (Hansen 1998) for details. In the inversions we use coarser meshes with \(256\times 128\) and \(96\times 48\times 48\) cells for the 2D instance and 3D instance, respectively. We consider the datasets generated in the preceding section and perturb the generated data with 10% iid Gaussian white noise.
To compute the L-curve, we solve 30 instances of the continuous relaxation for different values of \(\alpha\) that are logarithmically spaced between 1 and \(10^{-6}\). To accelerate the computation, we initialize the optimization with the solution from the previous \(\alpha\). For each value of \(\alpha\), we use up to 20 Gauss-Newton iterations and approximately compute the search direction using 5 iterations of projected preconditioned CG that use the Hessian of the regularization function as a preconditioner. For each experiment, we store the reconstructed source, the predicted data, the value of the misfit function, and the values of the regularization function (without the factor \(\alpha\)). The L-curve shown in Fig. 9 show the value of the regularizer and the value of the misfit of these optimal solutions. As is common, the axes are scaled logarithmically; and to provide additional insight, we have added visualizations of the reconstructed sources for the largest and smallest value of \(\alpha\) (resulting in overly smoothed and very noisy reconstructions, respectively) as well as solutions that provide a good trade-off. Using this process we select the regularization parameters \(\alpha =8.531 \times 10^{-3}\) for the two-dimensional instance and \(\alpha =5.298\times 10^{-3}\) for the three-dimensional instance, respectively. Computing the L-curves took about 4 and 48 minutes and involved 32,580 and 30,892 PDE solves in 2D and 3D, respectively. The large number of PDE solves underscores the importance of computing a factorization (or a good preconditioner in large-scale problems) apriori.
Rights and permissions
About this article
Cite this article
Sharma, M., Hahn, M., Leyffer, S. et al. Inversion of convection–diffusion equation with discrete sources. Optim Eng 22, 1419–1457 (2021). https://doi.org/10.1007/s11081-020-09536-5
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11081-020-09536-5
Keywords
- Mixed-integer optimization
- PDE-constrained optimization
- Mixed-integer nonlinear programming
- Source inversion