1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17 
18 
19 
20 #ifndef LIBMESH_OPTIMIZATION_SYSTEM_H
21 #define LIBMESH_OPTIMIZATION_SYSTEM_H
22 
23 // Local Includes
24 #include "libmesh/implicit_system.h"
25 
26 // C++ includes
27 
28 namespace libMesh
29 {
30 
31 
32 // Forward declarations
33 template<typename T> class OptimizationSolver;
34 
35 
36 /**
37  * This System subclass enables us to assemble an objective function,
38  * gradient, Hessian and bounds for optimization problems.
39  *
40  * \author David Knezevic
41  * \date 2015
42  */
43 class OptimizationSystem : public ImplicitSystem
44 {
45 public:
46 
47   /**
48    * Constructor.  Optionally initializes required
49    * data structures.
50    */
51   OptimizationSystem (EquationSystems & es,
52                       const std::string & name,
53                       const unsigned int number);
54 
55   /**
56    * Destructor.
57    */
58   virtual ~OptimizationSystem ();
59 
60   /**
61    * The type of system.
62    */
63   typedef OptimizationSystem sys_type;
64 
65   /**
66    * The type of the parent.
67    */
68   typedef ImplicitSystem Parent;
69 
70   /**
71    * Abstract base class to be used to calculate the objective
72    * function for optimization.
73    */
74   class ComputeObjective
75   {
76   public:
~ComputeObjective()77     virtual ~ComputeObjective () {}
78 
79     /**
80      * This function will be called to compute the objective function
81      * to be minimized, and must be implemented by the user in a
82      * derived class.
83      *
84      * \returns The value of the objective function at iterate \p X.
85      */
86     virtual Number objective (const NumericVector<Number> & X,
87                               sys_type & S) = 0;
88   };
89 
90 
91   /**
92    * Abstract base class to be used to calculate the gradient of
93    * an objective function.
94    */
95   class ComputeGradient
96   {
97   public:
~ComputeGradient()98     virtual ~ComputeGradient () {}
99 
100     /**
101      * This function will be called to compute the gradient of the
102      * objective function, and must be implemented by the user in
103      * a derived class. Set \p grad_f to be the gradient at the
104      * iterate \p X.
105      */
106     virtual void gradient (const NumericVector<Number> & X,
107                            NumericVector<Number> & grad_f,
108                            sys_type & S) = 0;
109   };
110 
111 
112   /**
113    * Abstract base class to be used to calculate the Hessian
114    * of an objective function.
115    */
116   class ComputeHessian
117   {
118   public:
~ComputeHessian()119     virtual ~ComputeHessian () {}
120 
121     /**
122      * This function will be called to compute the Hessian of
123      * the objective function, and must be implemented by the
124      * user in a derived class. Set \p H_f to be the gradient
125      * at the iterate \p X.
126      */
127     virtual void hessian (const NumericVector<Number> & X,
128                           SparseMatrix<Number> & H_f,
129                           sys_type & S) = 0;
130   };
131 
132   /**
133    * Abstract base class to be used to calculate the equality constraints.
134    */
135   class ComputeEqualityConstraints
136   {
137   public:
~ComputeEqualityConstraints()138     virtual ~ComputeEqualityConstraints () {}
139 
140     /**
141      * This function will be called to evaluate the equality constraints
142      * vector C_eq(X). This will impose the constraints C_eq(X) = 0.
143      */
144     virtual void equality_constraints (const NumericVector<Number> & X,
145                                        NumericVector<Number> & C_eq,
146                                        sys_type & S) = 0;
147   };
148 
149   /**
150    * Abstract base class to be used to calculate the Jacobian of the
151    * equality constraints.
152    */
153   class ComputeEqualityConstraintsJacobian
154   {
155   public:
~ComputeEqualityConstraintsJacobian()156     virtual ~ComputeEqualityConstraintsJacobian () {}
157 
158     /**
159      * This function will be called to evaluate the Jacobian of C_eq(X).
160      */
161     virtual void equality_constraints_jacobian (const NumericVector<Number> & X,
162                                                 SparseMatrix<Number> & C_eq_jac,
163                                                 sys_type & S) = 0;
164   };
165 
166   /**
167    * Abstract base class to be used to calculate the inequality constraints.
168    */
169   class ComputeInequalityConstraints
170   {
171   public:
~ComputeInequalityConstraints()172     virtual ~ComputeInequalityConstraints () {}
173 
174     /**
175      * This function will be called to evaluate the equality constraints
176      * vector C_ineq(X). This will impose the constraints C_ineq(X) >= 0.
177      */
178     virtual void inequality_constraints (const NumericVector<Number> & X,
179                                          NumericVector<Number> & C_ineq,
180                                          sys_type & S) = 0;
181   };
182 
183   /**
184    * Abstract base class to be used to calculate the Jacobian of the
185    * inequality constraints.
186    */
187   class ComputeInequalityConstraintsJacobian
188   {
189   public:
~ComputeInequalityConstraintsJacobian()190     virtual ~ComputeInequalityConstraintsJacobian () {}
191 
192     /**
193      * This function will be called to evaluate the Jacobian of C_ineq(X).
194      */
195     virtual void inequality_constraints_jacobian (const NumericVector<Number> & X,
196                                                   SparseMatrix<Number> & C_ineq_jac,
197                                                   sys_type & S) = 0;
198   };
199 
200   /**
201    * Abstract base class to be used to calculate the lower and upper
202    * bounds for all dofs in the system.
203    */
204   class ComputeLowerAndUpperBounds
205   {
206   public:
~ComputeLowerAndUpperBounds()207     virtual ~ComputeLowerAndUpperBounds () {}
208 
209     /**
210      * This function should update the following two vectors:
211      *   this->get_vector("lower_bounds"),
212      *   this->get_vector("upper_bounds").
213      */
214     virtual void lower_and_upper_bounds (sys_type & S) = 0;
215   };
216 
217   /**
218    * \returns A reference to *this.
219    */
system()220   sys_type & system () { return *this; }
221 
222   /**
223    * Clear all the data structures associated with
224    * the system.
225    */
226   virtual void clear () override;
227 
228   /**
229    * Initializes new data members of the system.
230    */
231   virtual void init_data () override;
232 
233   /**
234    * Reinitializes the member data fields associated with
235    * the system, so that, e.g., \p assemble() may be used.
236    */
237   virtual void reinit () override;
238 
239   /**
240    * Solves the optimization problem.
241    */
242   virtual void solve () override;
243 
244   /**
245    * Initialize storage for the equality constraints, and the
246    * corresponding Jacobian. The length of \p constraint_jac_sparsity
247    * indicates the number of constraints that will be imposed,
248    * and n_dofs_per_constraint[i] gives the indices that are non-zero
249    * in row i of the Jacobian.
250    */
251   void initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity);
252 
253   /**
254    * Initialize storage for the inequality constraints, as per
255    * initialize_equality_constraints_storage.
256    */
257   void initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity);
258 
259   /**
260    * \returns \p "Optimization".  Helps in identifying
261    * the system type in an equation system file.
262    */
system_type()263   virtual std::string system_type () const override { return "Optimization"; }
264 
265   /**
266    * The \p OptimizationSolver that is used for performing the optimization.
267    */
268   std::unique_ptr<OptimizationSolver<Number>> optimization_solver;
269 
270   /**
271    * The vector that stores equality constraints.
272    */
273   std::unique_ptr<NumericVector<Number>> C_eq;
274 
275   /**
276    * The sparse matrix that stores the Jacobian of C_eq.
277    */
278   std::unique_ptr<SparseMatrix<Number>> C_eq_jac;
279 
280   /**
281    * The vector that stores inequality constraints.
282    */
283   std::unique_ptr<NumericVector<Number>> C_ineq;
284 
285   /**
286    * The sparse matrix that stores the Jacobian of C_ineq.
287    */
288   std::unique_ptr<SparseMatrix<Number>> C_ineq_jac;
289 
290   /**
291    * Vectors to store the dual variables associated with equality
292    * and inequality constraints.
293    */
294   std::unique_ptr<NumericVector<Number>> lambda_eq;
295   std::unique_ptr<NumericVector<Number>> lambda_ineq;
296 
297   /**
298    * A copy of the equality and inequality constraint Jacobian sparsity
299    * patterns.
300    */
301   std::vector<std::set<numeric_index_type>> eq_constraint_jac_sparsity;
302   std::vector<std::set<numeric_index_type>> ineq_constraint_jac_sparsity;
303 };
304 
305 } // namespace libMesh
306 
307 #endif // LIBMESH_OPTIMIZATION_SYSTEM_H
308