MMA.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2020-2023 PCOpt/NTUA
    Copyright (C) 2020-2023 FOSS GP
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.


Class
    Foam::MMA

Description
    Update design variables using the Method of Moving Assymptotes.
    Can handle inequality constraints.

    Reference:
    \verbatim
        Svanberg, K. (1987)
        The method of moving asymptotes—a new method for structural
        optimization
        International Journal for Numerical Methods in Engineering, 24(2),
        359-373.
        https://doi.org/10.1002/nme.1620240207
    \endverbatim

SourceFiles
    MMA.C

\*---------------------------------------------------------------------------*/

#ifndef MMA_H
#define MMA_H

#include "constrainedOptimisationMethod.H"
#include "updateMethod.H"
#include "PtrList.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                             Class MMA Declaration
\*---------------------------------------------------------------------------*/

class MMA
:
    public constrainedOptimisationMethod,
    public updateMethod
{
protected:

    // Protected data

        //- The previous design variables. Used to compute new bounds
        scalarField x0_;

        //- The twice previous design variables. Used to compute new bounds
        scalarField x00_;

        //- The set of design variables being updated during the subproblem.
        //  Correction will end up being the difference of this and the
        //  old values of the design variables
        scalarField xNew_;

        //- Old objective value
        //  Needed by GCMMA
        scalar oldObjectiveValue_;

        //- Old constraint values
        //  Needed by GCMMA
        scalarField oldCValues_;

        //- Objective/Constraint values and approximations in the new point
        //  Needed by GCMMA
        PtrList<scalarField> valsAndApproxs_;

        //- Second term in the approximation function
        scalar z_;

        //- Term multiplying z in the objective function
        scalar alpha0_;

        //- Terms multyplying z in the constraint functions
        scalarField alpha_;

        //- y in the constraints functions
        scalarField y_;

        //- y multipliers in the objective function
        scalarField c_;

        //- y^2 multipliers in the objective function
        scalarField d_;

        //- Lower design variables asymptotes
        scalarField lower_;

        //- Upper design variables asymptotes
        scalarField upper_;

        //- Lower design variables constraints
        scalarField a_;

        //- Upper design variables constraints
        scalarField b_;

        //- Constants in the (p,q) functions.
        //  Size is cValues_.size() + 1 (the objective)
        scalarField rho_;

        //- Bound the rho value with an upper value?
        bool boundRho_;

        //- Correct the design variables
        bool correctDVs_;

        // Dual variables

            //- Constraint Lagrange multipliers
            scalarField lamda_;

            //- Lagrange multipliers for the lower limits constraints
            scalarField ksi_;

            //- Lagrange multipliers for the upper limits constraints
            scalarField Eta_;

            //- Lagrange multipliers for the y constraints
            scalarField mu_;

            //- Lagrange multiplier for the z constraint
            scalar zeta_;

            //- Slack variables for the inequality constraints
            scalarField s_;


        // Fields holding updates of the design, Lagrange and slack variables

            scalarField deltaLamda_;
            scalar deltaZ_;
            scalarField deltaX_;
            scalarField deltaY_;
            scalarField deltaS_;
            scalar deltaZeta_;
            scalarField deltaMu_;
            scalarField deltaEta_;
            scalarField deltaKsi_;


        //- Infinitesimal quantity
        //  Updated during the inner iterations of the subproblem
        scalar eps_;

        //- Maxmimum number of Newton iterations for the subproblem
        label maxNewtonIters_;

        //- Maxmimum number of line search iterations for each iteration of the
        //- subproblem
        label maxLineSearchIters_;

        //- Initialize every subproblem with x = 0.5*(a + b)
        //- and reset Lagrange multipliers
        //  The solution of the subproblem can become relatively expensive
        //  in some cases, hence the GCMMA subproblems are usually
        //  initialized with the solution of the previous inner iteration,
        //  to save up some time. Occassionally, this can cause the solution
        //  of the subproblem to be trapped in a non-feasible point.
        //  In such cases, the solution of each subproblem could be
        //  initialized from scratch
        bool initializeEverySubproblem_;


        // Constant parameters

            //- Lower assymprote reduction rate
            scalar asymDecr_;

            //- Upper assymprote increase rate
            scalar asymIncr_;

            //- Used to update the assymptotes in the first optimisation cycle
            scalar sInit_;

            //- Movement of the unatainable upper and lower bounds
            scalar move_;

            //- Constant in p, q functions
            scalar raa0_;


        // Parameters related to the update of rho in GCMMA

            //- Multiplier of the mean derivative used for the initialisation
            //- of rho
            scalar maxInitRhoMult_;

            //- Multiplier for the max. rho value during its update
            scalar maxRhoMult_;

            //- Change rho constants in each iteration?
            bool variableRho_;

            //- Change rho constants in each iteration?
            bool dynamicRhoInitialisation_;

            //- The rho values obtained by the dynamic rho initialisation
            //- might be too conservative. Multiply with this number to relax
            //- them
            scalar dynamicRhoMult_;


        // Normalisation based on the gradient

            //- Perform the normalisation
            bool normalise_;

            //- Perform the normalisation in each optimisation cycle
            bool continuousNormalisation_;

            //- Normalisation factor for the objective
            autoPtr<scalar> Jnorm_;

            //- Normalisation factors for the constraints
            tmp<scalarField> Cnorm_;

            //- Constaint weight after the normalisation
            scalar cw_;

            //- Constaint weight after the normalisation
            label lastNormalisationStep_;


    // Protected Member Functions

        //- Update sizes of fields related to the active design variables
        void updateSizes();

        //- Initialize rho constants in the (p, q) functions
        //  These control how aggressive or conservative the method will be
        void initializeRho();

        //- Update the bounds associated with the design variables
        void updateBounds();

        //- Allocate fields related to constraints
        void initialize();

        //- Store old objcective and constraint values
        //  Needed by GCMMA
        void storeOldValues();

        //- p and q functions, used to approximate the objective and contraints
        //  Computed based on the current set of design variables, not updated
        //  during the iterations of the subproblem
        tmp<scalarField> p
        (
            const scalarField& derivs,
            const scalar r,
            const scalarField& x
        );
        tmp<scalarField> p(const scalarField& derivs, const scalar r);
        tmp<scalarField> q
        (
            const scalarField& derivs,
            const scalar r,
            const scalarField& x
        );
        tmp<scalarField> q(const scalarField& derivs, const scalar r);

        //- g of all constraint functions
        //  Computed using the current set of design variables, updated during
        //  the solution of the subproblem
        tmp<scalarField> gConstr(const scalarField& vars);

        //- The rhs of the contraint approximations.
        //  Always computed using the current (frozen) set of design variables
        tmp<scalarField> b();

        //- p and q functions, using the dual lamda
        tmp<scalarField> pLamda();
        tmp<scalarField> qLamda();

        //- Zero the updates computed in the previous optimisation cycle
        void zeroUpdates();

        //- Computation of the approximate function
        scalar fTilda
        (
            const scalarField& xInit,
            const scalarField& x,
            const scalar f,
            const scalarField& dfdx,
            const scalar rho
        );

        //- Read in scalarField from self (potentially in binary), or allocate
        //- field with the size of the active design variables and given value
        tmp<scalarField> getOrDefaultScalarField
        (
            const word& name,
            const label size,
            const scalar value = Zero
        );

        //- Read in scalarField from self (potentially in binary), or allocate
        //- field with the size of the active design variables and given value
        void setOrDefaultScalarField
        (
            scalarField& field,
            const word& name,
            const label size,
            const scalarField& defaultField
        );


        // Functions related to the solution of the primal-dual subproblem

            //- Compute direction for the Newton problem
            void computeNewtonDirection();

            //- Perform line search and return max residual corresponding to
            //- the updated solution
            scalar lineSearch();

            //- Update the current solution using the known direction and the
            //- given step length
            void updateSolution(const scalar step);

            //- Adjust step to satisfy cireteria
            void adjustStep
            (
                scalar& step,
                const scalar value,
                const scalar update
            );

            //- Compute and return residuals based on the current solution
            tmp<scalarField> computeResiduals();


        //- Normalise the objective and constraints if necessary
        void normalise();


private:

    // Private Member Functions

        //- No copy construct
        MMA(const MMA&);

        //- No copy assignment
        void operator=(const MMA&);


public:

    //- Runtime type information
    TypeName("MMA");


    // Constructors

        //- Construct from components
        MMA
        (
            const fvMesh& mesh,
            const dictionary& dict,
            autoPtr<designVariables>& designVars,
            const label nConstraints,
            const word& type
        );


    //- Destructor
    virtual ~MMA() = default;


    // Member Functions

        //- Compute design variables correction
        void computeCorrection();

        //- Solve subproblem using a Newton optimiser
        void solveSubproblem();

        // Functions needed by GCMMA

            //- Compute objective/constraint values and their approximations
            //  Needed to check convergence and update rho, if necessary
            void updateValuesAndApproximations();

            //- Get objective/constraint values and their approximations
            const PtrList<scalarField>& getValuesAndApproximations() const;

            //- Update rho value if needed
            void updateRho();

            //- Get rho value if needed
            const scalarField& getRho() const;

            //- Set variable rho
            void setVariableRho(bool varRho = true);

            //- Checks convergence of GCMMA
            bool converged();


        //- Write useful quantities to files
        virtual bool writeData(Ostream& os) const;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //