47 constrainedOptimisationMethod,
56 void Foam::SQP::updateHessian()
76 Info<<
"Scaling Hessian with factor " << scaleFactor <<
endl;
79 Hessian_()[varI][varI] /= scaleFactor;
85 <<
" y*s is negative. Skipping the scaling of the first Hessian" 96 <<
" y*s is below threshold. Using damped form" <<
endl;
101 <<
"Unmodified Hessian curvature index " << ys <<
endl;
112 void Foam::SQP::update()
115 SquareMatrix<scalar> HessianInv =
inv(Hessian_());
118 Info<<
"Hessian " << Hessian_() <<
endl;
119 Info<<
"HessianInv " << HessianInv <<
endl;
120 label
n = Hessian_().n();
121 SquareMatrix<scalar> test(
n,
Zero);
122 for (label
k = 0;
k <
n;
k++)
124 for (label l = 0; l <
n; l++)
127 for (label i = 0; i <
n; i++)
129 elem += Hessian_()[
k][i] * HessianInv[i][l];
134 Info<<
"Validation " << test <<
endl;
138 label nc = constraintDerivatives_.size();
142 activeDerivs.
map(LagrangianDerivatives_, activeDesignVars_);
143 scalarField WgradL = rightMult(HessianInv, activeDerivs);
149 activeConsDerivs.
map(constraintDerivatives_[cI], activeDesignVars_);
150 lamdaRHS[cI] = globalSum(activeConsDerivs * WgradL) - cValues_[cI];
153 Info<<
"lamdaRHS total|deriv part|constraint part " 154 << lamdaRHS[cI] <<
" " << globalSum(activeConsDerivs * WgradL)
155 <<
" " << cValues_[cI] <<
endl;
160 SquareMatrix<scalar> AWA(nc,
Zero);
161 PtrList<scalarField> WA(nc);
162 for (label j = 0; j < nc; j++)
165 gradcJ.
map(constraintDerivatives_[j], activeDesignVars_);
166 WA.set(j,
new scalarField(rightMult(HessianInv, gradcJ)));
167 for (label i = 0; i < nc; i++)
170 gradcI.
map(constraintDerivatives_[i], activeDesignVars_);
171 AWA[i][j] = globalSum(gradcI * WA[j]);
174 SquareMatrix<scalar> invAWA =
inv(AWA);
175 scalarField deltaLamda = rightMult(invAWA, lamdaRHS);
180 Info<<
"lamda update " << deltaLamda <<
endl;
182 lamdas_ += deltaLamda;
188 activeCorrection += WA[cI]*deltaLamda[cI];
190 activeCorrection *= etaHessian_;
193 forAll(activeDesignVars_, varI)
195 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
204 void Foam::SQP::storeOldFields()
206 derivativesOld_ = objectiveDerivatives_;
207 if (constraintDerivativesOld_.empty())
209 constraintDerivativesOld_.setSize(constraintDerivatives_.size());
211 forAll(constraintDerivativesOld_, cI)
213 constraintDerivativesOld_[cI] = constraintDerivatives_[cI];
215 correctionOld_ = correction_;
219 Foam::scalar Foam::SQP::meritFunctionConstraintPart()
const 236 const label nConstraints,
244 coeffsDict(
type).getOrDefault<scalar>(
"dumpingThreshold", 0.2)
255 LagrangianDerivatives_ = objectiveDerivatives_;
266 if (mu_ <
max(
mag(lamdas_)) + delta_)
268 mu_ =
max(
mag(lamdas_)) + 2*delta_;
271 Info<<
"Updated mu value to " << mu_ <<
endl;
274 scalar
L = objectiveValue_ + mu_*
sum(
mag(cValues_));
283 globalSum(objectiveDerivatives_*correction_)
void size(const label n)
Older name for setAddressableSize.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A list of keyword definitions, which are a keyword followed by a number of values (eg...
const vector L(dict.get< vector >("L"))
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual bool addToFile(Ostream &os) const
Write continuation info.
bool scaleFirstHessian_
Scale the initial unitary Hessian approximation.
label k
Boltzmann constant.
virtual bool writeData(Ostream &os) const
Write useful quantities to files.
const scalarField rightMult(const SquareMatrix< scalar > &, const scalarField &)
const scalarField leftMult(const scalarField &, const SquareMatrix< scalar > &)
scalar dumpingThreshold_
Curvature threshold.
scalarField LagrangianDerivatives_
Derivatives of the Lagrangian function.
Macros for easy insertion into run-time selection tables.
PtrList< scalarField > constraintDerivatives_
Derivatives of the constraints.
#define forAll(list, i)
Loop across all elements in list.
void computeCorrection()
Compute design variables correction.
void computeCorrection()
Compute design variables correction.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const labelList & activeDesignVars_
Map to active design variables.
Base class for quasi-Newton methods.
A class for handling words, derived from Foam::string.
autoPtr< SquareMatrix< scalar > > Hessian_
The Hessian or its inverse, depending on the deriving class.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label counter_
Optimisation cycle count.
scalarField derivativesOld_
The previous derivatives.
List< scalarField > constraintDerivativesOld_
The previous constraint derivatives.
virtual bool writeAuxiliaryData()
Write merit function information.
Base class for Sequantial Quadratic Programming (SQP) methods.
#define DebugInfo
Report an information message using Foam::Info.
void map(const UList< Type > &mapF, const labelUList &mapAddressing)
1 to 1 map from the given field
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Istream and Ostream manipulators taking arguments.
scalar globalSum(const scalarField &field)
Compute either global or local sum, based on globalSum flag.
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
virtual bool writeMeritFunction(const updateMethod &UpdateMethod)
Write info about the merit function.
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh data needed to do the Finite Volume discretisation.
SquareMatrix< scalar > outerProd(const scalarField &, const scalarField &)
virtual scalar meritFunctionDirectionalDerivative()
Derivative of the merit function. Could be different than the objective derivative in the presence of...
messageStream Info
Information stream (stdout output on master, null elsewhere)
scalarField lamdas_
Lagrange multipliers.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
virtual scalar computeMeritFunction()
Compute merit function. Could be different than the objective in the presence of constraints.
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
virtual bool writeData(Ostream &os) const
Write useful quantities to files.
static constexpr const zero Zero
Global zero (0)
scalarField correctionOld_
The previous correction.