Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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          10 : class OptimizationSystem : public ImplicitSystem
      44             : {
      45             : public:
      46             : 
      47             :   /**
      48             :    * Constructor.
      49             :    */
      50             :   OptimizationSystem (EquationSystems & es,
      51             :                       const std::string & name,
      52             :                       const unsigned int number);
      53             : 
      54             :   /**
      55             :    * Special functions.
      56             :    * - This class has the same restrictions/defaults as its base class.
      57             :    * - The destructor is defaulted out-of-line.
      58             :    */
      59             :   OptimizationSystem (const OptimizationSystem &) = delete;
      60             :   OptimizationSystem & operator= (const OptimizationSystem &) = delete;
      61             :   OptimizationSystem (OptimizationSystem &&) = default;
      62             :   OptimizationSystem & operator= (OptimizationSystem &&) = delete;
      63             :   virtual ~OptimizationSystem ();
      64             : 
      65             :   /**
      66             :    * The type of system.
      67             :    */
      68             :   typedef OptimizationSystem sys_type;
      69             : 
      70             :   /**
      71             :    * The type of the parent.
      72             :    */
      73             :   typedef ImplicitSystem Parent;
      74             : 
      75             :   /**
      76             :    * Abstract base class to be used to calculate the objective
      77             :    * function for optimization.
      78             :    */
      79             :   class ComputeObjective
      80             :   {
      81             :   public:
      82             :     virtual ~ComputeObjective () = default;
      83             : 
      84             :     /**
      85             :      * This function will be called to compute the objective function
      86             :      * to be minimized, and must be implemented by the user in a
      87             :      * derived class.
      88             :      *
      89             :      * \returns The value of the objective function at iterate \p X.
      90             :      */
      91             :     virtual Number objective (const NumericVector<Number> & X,
      92             :                               sys_type & S) = 0;
      93             :   };
      94             : 
      95             : 
      96             :   /**
      97             :    * Abstract base class to be used to calculate the gradient of
      98             :    * an objective function.
      99             :    */
     100             :   class ComputeGradient
     101             :   {
     102             :   public:
     103             :     virtual ~ComputeGradient () = default;
     104             : 
     105             :     /**
     106             :      * This function will be called to compute the gradient of the
     107             :      * objective function, and must be implemented by the user in
     108             :      * a derived class. Set \p grad_f to be the gradient at the
     109             :      * iterate \p X.
     110             :      */
     111             :     virtual void gradient (const NumericVector<Number> & X,
     112             :                            NumericVector<Number> & grad_f,
     113             :                            sys_type & S) = 0;
     114             :   };
     115             : 
     116             : 
     117             :   /**
     118             :    * Abstract base class to be used to calculate the Hessian
     119             :    * of an objective function.
     120             :    */
     121             :   class ComputeHessian
     122             :   {
     123             :   public:
     124             :     virtual ~ComputeHessian () = default;
     125             : 
     126             :     /**
     127             :      * This function will be called to compute the Hessian of
     128             :      * the objective function, and must be implemented by the
     129             :      * user in a derived class. Set \p H_f to be the gradient
     130             :      * at the iterate \p X.
     131             :      */
     132             :     virtual void hessian (const NumericVector<Number> & X,
     133             :                           SparseMatrix<Number> & H_f,
     134             :                           sys_type & S) = 0;
     135             :   };
     136             : 
     137             :   /**
     138             :    * Abstract base class to be used to calculate the equality constraints.
     139             :    */
     140             :   class ComputeEqualityConstraints
     141             :   {
     142             :   public:
     143             :     virtual ~ComputeEqualityConstraints () = default;
     144             : 
     145             :     /**
     146             :      * This function will be called to evaluate the equality constraints
     147             :      * vector C_eq(X). This will impose the constraints C_eq(X) = 0.
     148             :      */
     149             :     virtual void equality_constraints (const NumericVector<Number> & X,
     150             :                                        NumericVector<Number> & C_eq,
     151             :                                        sys_type & S) = 0;
     152             :   };
     153             : 
     154             :   /**
     155             :    * Abstract base class to be used to calculate the Jacobian of the
     156             :    * equality constraints.
     157             :    */
     158             :   class ComputeEqualityConstraintsJacobian
     159             :   {
     160             :   public:
     161             :     virtual ~ComputeEqualityConstraintsJacobian () = default;
     162             : 
     163             :     /**
     164             :      * This function will be called to evaluate the Jacobian of C_eq(X).
     165             :      */
     166             :     virtual void equality_constraints_jacobian (const NumericVector<Number> & X,
     167             :                                                 SparseMatrix<Number> & C_eq_jac,
     168             :                                                 sys_type & S) = 0;
     169             :   };
     170             : 
     171             :   /**
     172             :    * Abstract base class to be used to calculate the inequality constraints.
     173             :    */
     174             :   class ComputeInequalityConstraints
     175             :   {
     176             :   public:
     177             :     virtual ~ComputeInequalityConstraints () = default;
     178             : 
     179             :     /**
     180             :      * This function will be called to evaluate the equality constraints
     181             :      * vector C_ineq(X). This will impose the constraints C_ineq(X) >= 0.
     182             :      */
     183             :     virtual void inequality_constraints (const NumericVector<Number> & X,
     184             :                                          NumericVector<Number> & C_ineq,
     185             :                                          sys_type & S) = 0;
     186             :   };
     187             : 
     188             :   /**
     189             :    * Abstract base class to be used to calculate the Jacobian of the
     190             :    * inequality constraints.
     191             :    */
     192             :   class ComputeInequalityConstraintsJacobian
     193             :   {
     194             :   public:
     195             :     virtual ~ComputeInequalityConstraintsJacobian () = default;
     196             : 
     197             :     /**
     198             :      * This function will be called to evaluate the Jacobian of C_ineq(X).
     199             :      */
     200             :     virtual void inequality_constraints_jacobian (const NumericVector<Number> & X,
     201             :                                                   SparseMatrix<Number> & C_ineq_jac,
     202             :                                                   sys_type & S) = 0;
     203             :   };
     204             : 
     205             :   /**
     206             :    * Abstract base class to be used to calculate the lower and upper
     207             :    * bounds for all dofs in the system.
     208             :    */
     209             :   class ComputeLowerAndUpperBounds
     210             :   {
     211             :   public:
     212             :     virtual ~ComputeLowerAndUpperBounds () = default;
     213             : 
     214             :     /**
     215             :      * This function should update the following two vectors:
     216             :      *   this->get_vector("lower_bounds"),
     217             :      *   this->get_vector("upper_bounds").
     218             :      */
     219             :     virtual void lower_and_upper_bounds (sys_type & S) = 0;
     220             :   };
     221             : 
     222             :   /**
     223             :    * \returns A reference to *this.
     224             :    */
     225             :   sys_type & system () { return *this; }
     226             : 
     227             :   /**
     228             :    * Clear all the data structures associated with
     229             :    * the system.
     230             :    */
     231             :   virtual void clear () override;
     232             : 
     233             :   /**
     234             :    * Initializes new data members of the system.
     235             :    */
     236             :   virtual void init_data () override;
     237             : 
     238             :   /**
     239             :    * Reinitializes the member data fields associated with
     240             :    * the system, so that, e.g., \p assemble() may be used.
     241             :    */
     242             :   virtual void reinit () override;
     243             : 
     244             :   /**
     245             :    * Solves the optimization problem.
     246             :    */
     247             :   virtual void solve () override;
     248             : 
     249             :   /**
     250             :    * Initialize storage for the equality constraints, and the
     251             :    * corresponding Jacobian. The length of \p constraint_jac_sparsity
     252             :    * indicates the number of constraints that will be imposed,
     253             :    * and n_dofs_per_constraint[i] gives the indices that are non-zero
     254             :    * in row i of the Jacobian.
     255             :    */
     256             :   void initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity);
     257             : 
     258             :   /**
     259             :    * Initialize storage for the inequality constraints, as per
     260             :    * initialize_equality_constraints_storage.
     261             :    */
     262             :   void initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity);
     263             : 
     264             :   /**
     265             :    * \returns \p "Optimization".  Helps in identifying
     266             :    * the system type in an equation system file.
     267             :    */
     268          71 :   virtual std::string system_type () const override { return "Optimization"; }
     269             : 
     270             :   /**
     271             :    * The \p OptimizationSolver that is used for performing the optimization.
     272             :    */
     273             :   std::unique_ptr<OptimizationSolver<Number>> optimization_solver;
     274             : 
     275             :   /**
     276             :    * The vector that stores equality constraints.
     277             :    */
     278             :   std::unique_ptr<NumericVector<Number>> C_eq;
     279             : 
     280             :   /**
     281             :    * The sparse matrix that stores the Jacobian of C_eq.
     282             :    */
     283             :   std::unique_ptr<SparseMatrix<Number>> C_eq_jac;
     284             : 
     285             :   /**
     286             :    * The vector that stores inequality constraints.
     287             :    */
     288             :   std::unique_ptr<NumericVector<Number>> C_ineq;
     289             : 
     290             :   /**
     291             :    * The sparse matrix that stores the Jacobian of C_ineq.
     292             :    */
     293             :   std::unique_ptr<SparseMatrix<Number>> C_ineq_jac;
     294             : 
     295             :   /**
     296             :    * Vectors to store the dual variables associated with equality
     297             :    * and inequality constraints.
     298             :    */
     299             :   std::unique_ptr<NumericVector<Number>> lambda_eq;
     300             :   std::unique_ptr<NumericVector<Number>> lambda_ineq;
     301             : 
     302             :   /**
     303             :    * A copy of the equality and inequality constraint Jacobian sparsity
     304             :    * patterns.
     305             :    */
     306             :   std::vector<std::set<numeric_index_type>> eq_constraint_jac_sparsity;
     307             :   std::vector<std::set<numeric_index_type>> ineq_constraint_jac_sparsity;
     308             : };
     309             : 
     310             : } // namespace libMesh
     311             : 
     312             : #endif // LIBMESH_OPTIMIZATION_SYSTEM_H