46 constrainedOptimisationMethod,
57 { preconditioners::invHessian,
"invHessian" },
58 { preconditioners::ShermanMorrison,
"ShermanMorrison" }
99 if (includeBoundConstraints_)
102 const label
n(activeDesignVars_.size());
126 const label m(nConstraints_);
128 if (includeExtraVars_)
137 doAllocateLagrangeMultipliers_ =
false;
144 scalarField LagrangianDerivativesOld(derivativesOld_);
145 forAll(constraintDerivatives_, cI)
147 LagrangianDerivatives_ += lamdas_[cI]*constraintDerivatives_[cI];
148 LagrangianDerivativesOld += lamdas_[cI]*constraintDerivativesOld_[cI];
151 if (includeBoundConstraints_)
153 forAll(activeDesignVars_, aI)
155 const label varI(activeDesignVars_[aI]);
156 const scalar contr(uTilda_()[aI] - lTilda_()[aI]);
157 LagrangianDerivatives_[varI] += contr;
158 LagrangianDerivativesOld[varI] += contr;
163 updateVectors(LagrangianDerivatives_, LagrangianDerivativesOld);
169 const scalarField x(designVars_().getVars(), activeDesignVars_);
173 if (includeBoundConstraints_)
175 lTilda_() = scalar(1);
176 uTilda_() = scalar(1);
185 if (includeExtraVars_)
187 extraVars_() = scalar(1);
188 z_() =
max(1, 0.5*c_);
202 if (includeBoundConstraints_)
204 deltaLTilda_() =
Zero;
206 deltaUTilda_() =
Zero;
210 if (includeExtraVars_)
212 deltaExtraVars_() =
Zero;
223 if (includeBoundConstraints_)
226 (uTilda_()*resFus() + resFuTilda())/us_()
227 - (lTilda_()*resFls() + resFlTilda())/ls_();
231 if (includeExtraVars_)
233 mult += extraVars_()/z_();
234 AMult -= (extraVars_()*resFExtraVars() + resFz())/z_();
239 const label varI(activeDesignVars_[aI]);
240 forAll(constraintDerivatives_, cI)
242 FDx[aI] += constraintDerivatives_[cI][varI]*AMult[cI];
259 forAll(constraintDerivatives_, cI)
261 const scalarField& cDerivsI = constraintDerivatives_[cI];
263 globalSum(
scalarField(cDerivsI, activeDesignVars_)*deltaP_);
265 GsLADp *= lamdas_/gs_;
270 const label varI(activeDesignVars_[aI]);
271 forAll(constraintDerivatives_, cI)
273 rhs[aI] += constraintDerivatives_[cI][varI]*GsLADp[cI];
278 if (includeBoundConstraints_)
280 rhs += (lTilda_()/ls_() + uTilda_()/us_())*deltaP_;
283 rhs = -invHessianVectorProduct(rhs);
285 rhs = 0.95*deltaP_ + 0.05*rhs;
297 scalar res(
sqrt(globalSum(r*r)));
299 scalar rz(globalSum(r*z));
305 scalar a = rz/globalSum(
p*Ap);
308 res =
sqrt(globalSum(r*r));
309 z = preconVectorProduct(r, precon);
311 scalar
beta = rz/rzOld;
314 }
while (iter++ < maxDxIters_ && res > relTol_*resInit);
316 <<
"CG, Solving for deltaP, Initial Residual " << resInit
317 <<
", Final Residual " << res
318 <<
", No Iterations " << iter <<
endl;
331 forAll(constraintDerivatives_, cI)
333 const scalarField& cDerivsI = constraintDerivatives_[cI];
338 if (includeExtraVars_)
340 mult += extraVars_()/z_();
347 const label varI(activeDesignVars_[aI]);
348 forAll(constraintDerivatives_, cI)
350 Ap[aI] += constraintDerivatives_[cI][varI]*GsLAv[cI];
355 if (includeBoundConstraints_)
357 Ap += (lTilda_()/ls_() + uTilda_()/us_())*
vector;
376 precon.
reset(diagPreconditioner().ptr());
380 else if (preconType_ == preconditioners::invHessian)
382 return invHessianVectorProduct(
vector);
384 else if (preconType_ == preconditioners::ShermanMorrison)
386 return ShermanMorrisonPrecon(
vector);
395 tmp<scalarField> tpreconditioner(HessianDiag());
396 scalarField& preconditioner = tpreconditioner.ref();
399 forAll(constraintDerivatives_, cI)
401 scalarField cDerivs(constraintDerivatives_[cI], activeDesignVars_);
402 scalar mult(gs_[cI]/lamdas_[cI]);
403 if (includeExtraVars_)
405 mult += extraVars_()[cI]/z_()[cI];
407 preconditioner +=
sqr(cDerivs)/mult;
410 if (includeBoundConstraints_)
412 preconditioner += lTilda_()/ls_() + uTilda_()/us_();
415 preconditioner = 1./preconditioner;
417 return tpreconditioner;
435 if (includeBoundConstraints_)
437 tdiag.
reset((lTilda_()/ls_() + uTilda_()/us_()).ptr());
442 PtrList<scalarField> r1Updates(cValues_.size());
444 forAll(constraintDerivatives_, cI)
446 const scalarField& cDerivsI = constraintDerivatives_[cI];
451 if (includeExtraVars_)
453 mult += extraVars_()/z_();
457 ShermanMorrisonRank1Update(r1Updates,
vector, tdiag, mult, mult.
size());
475 Ap = invHessianVectorProduct(
p, counter_,
diag.shallowClone());
482 Ap = ShermanMorrisonRank1Update(r1Updates,
p,
diag, mult,
n);
483 tmp<scalarField> tAv =
484 ShermanMorrisonRank1Update(r1Updates, r1Updates[
n],
diag, mult,
n);
486 scalar yHs = globalSum(r1Updates[
n]*Av)/mult[
n];
487 Ap -= Av*globalSum(r1Updates[
n]*Ap)/(1 + yHs)/mult[
n];
504 forAll(constraintDerivatives_, cI)
506 const scalarField& cDerivsI = constraintDerivatives_[cI];
508 globalSum(
scalarField(cDerivsI, activeDesignVars_)*deltaP_);
511 if (includeExtraVars_)
513 mult += extraVars_()/z_();
514 deltaLamda_ += (resFz() + extraVars_()*resFExtraVars())/z_();
516 deltaLamda_ += resFGs() - resFlamda()/lamdas_;
520 deltaGs_ = -(gs_*deltaLamda_ + resFlamda())/lamdas_;
522 if (includeBoundConstraints_)
525 deltaLs_() = deltaP_ + resFls();
528 deltaUs_() = -deltaP_ + resFus();
531 deltaLTilda_() = -(lTilda_()*deltaLs_() + resFlTilda())/ls_();
534 deltaUTilda_() = -(uTilda_()*deltaUs_() + resFuTilda())/us_();
537 if (includeExtraVars_)
539 deltaZ_() = -deltaLamda_ + resFExtraVars();
540 deltaExtraVars_() = - (extraVars_()*deltaZ_() + resFz())/z_();
547 const label
n(p_.size());
548 const label m(cValues_.size());
551 if (includeBoundConstraints_)
553 for (label i = 0; i <
n; ++i)
555 adjustStep(step, ls_()[i], deltaLs_()[i]);
556 adjustStep(step, us_()[i], deltaUs_()[i]);
557 adjustStep(step, lTilda_()[i], deltaLTilda_()[i]);
558 adjustStep(step, uTilda_()[i], deltaUTilda_()[i]);
563 for (label i = 0; i < m; ++i)
565 adjustStep(step, lamdas_[i], deltaLamda_[i]);
566 adjustStep(step, gs_[i], deltaGs_[i]);
567 if (includeExtraVars_)
569 adjustStep(step, extraVars_()[i], deltaExtraVars_()[i]);
570 adjustStep(step, z_()[i], deltaZ_()[i]);
578 reduce(step, minOp<scalar>());
585 Info<<
"Step before line search is " << step <<
endl;
589 scalar normResOld =
sqrt(globalSum(
magSqr(computeResiduals())));
590 scalar maxRes(GREAT);
592 for (label i = 0; i < maxLineSearchIters_ ; ++i)
595 updateSolution(step);
599 scalar normResNew =
sqrt(globalSum(
magSqr(resNew)));
602 if (normResNew < normResOld)
605 <<
"Initial residual = " << normResOld <<
", " 606 <<
"Final residual = " << normResNew <<
", " 607 <<
"No of LineSearch Iterations = " << i + 1
614 if (i != maxLineSearchIters_ - 1)
616 updateSolution(-step);
621 Info<<
tab <<
"Line search did not converge. " 622 <<
"Procceding with the last compute step" <<
endl;
629 Info<<
"Step after line search is " << step <<
nl <<
endl;
643 if (0.99*value + step*
update < scalar(0))
645 step = -0.99*value/
update;
653 lamdas_ += step*deltaLamda_;
654 gs_ += step*deltaGs_;
655 if (includeBoundConstraints_)
657 lTilda_() += step*deltaLTilda_();
658 ls_() += step*deltaLs_();
659 uTilda_() += step*deltaUTilda_();
660 us_() += step*deltaUs_();
662 if (includeExtraVars_)
664 extraVars_() += step*deltaExtraVars_();
665 z_() += step*deltaZ_();
672 const label
n(activeDesignVars_.size());
673 const label m(cValues_.size());
674 label size(includeBoundConstraints_ ? 5*
n + 2*m :
n + 2*m);
675 if (includeExtraVars_)
689 res.rmap(resFGs()(),
identity(m, iRes));
693 res.rmap(resFlamda()(),
identity(m, iRes));
696 if (includeBoundConstraints_)
707 res.rmap(resFlTilda()(),
identity(
n, iRes));
711 res.rmap(resFuTilda()(),
identity(
n, iRes));
715 if (includeExtraVars_)
718 res.rmap(resFExtraVars()(),
identity(m, iRes));
722 res.rmap(resFz(),
identity(m, iRes));
732 tmp<scalarField> tgradL
746 forAll(constraintDerivatives_, cI)
750 *
scalarField(constraintDerivatives_[cI], activeDesignVars_);
753 if (includeBoundConstraints_)
755 gradL += uTilda_() - lTilda_();
764 tmp<scalarField> tinvHFL
768 forAll(constraintDerivatives_, cI)
772 *
scalarField(constraintDerivatives_[cI], activeDesignVars_);
775 if (includeBoundConstraints_)
777 invHFL += uTilda_() - lTilda_();
780 invHFL = invHessianVectorProduct(invHFL);
795 forAll(constraintDerivatives_, cI)
800 scalarField(constraintDerivatives_[cI], activeDesignVars_)*p_
804 if (includeExtraVars_)
815 return (lamdas_*gs_ - eps_);
821 if (includeBoundConstraints_)
823 const scalarField x(designVars_().getVars(), activeDesignVars_);
825 (designVars_().lowerBounds()(), activeDesignVars_);
827 return (
x + p_ -
xMin - ls_());
835 if (includeBoundConstraints_)
837 const scalarField x(designVars_().getVars(), activeDesignVars_);
839 (designVars_().upperBounds()(), activeDesignVars_);
841 return (
xMax -
x - p_ - us_());
849 if (includeBoundConstraints_)
851 return (lTilda_()*ls_() - eps_);
859 if (includeBoundConstraints_)
861 return (uTilda_()*us_() - eps_);
869 if (includeExtraVars_)
871 return (c_ - lamdas_ - z_());
879 if (includeExtraVars_)
881 return (z_()*extraVars_() - eps_);
891 if (includeBoundConstraints_ || !cValues_.empty())
893 scalar resMax(
gMax(
mag(computeResiduals())));
898 <<
"Newton iter " << iter <<
nl <<
endl;
901 if (resMax < 0.9*eps_)
907 computeNewtonDirection();
912 resMax = lineSearch();
914 <<
"max residual = " << resMax <<
", " 915 <<
"eps = " << eps_ <<
nl <<
endl;
918 iter++ < maxNewtonIters_ && (eps_ > epsMin_ || resMax > 0.9*eps_)
920 Info<<
"Finished solving the QP subproblem" <<
nl <<
endl;
921 if (iter == maxNewtonIters_)
924 <<
"Iterative solution of the QP problem did not converge" 929 scalarField vars(designVars_().getVars(), activeDesignVars_);
938 if (includeBoundConstraints_)
944 Info<<
"Max of lTilda*ls " <<
gMax(lTilda_()*ls_()) <<
endl;
945 Info<<
"Max of uTilda*us " <<
gMax(uTilda_()*us_()) <<
endl;
947 if (includeExtraVars_)
949 Info<<
"Min of extraVars " <<
gMin(extraVars_()) <<
endl;
950 Info<<
"Max of extraVars*z " <<
gMax(extraVars_()*z_()) <<
endl;
956 computeNewtonDirection();
961 correction_.rmap(p_, activeDesignVars_);
968 correction_ *= etaHessian_;
975 derivativesOld_ = objectiveDerivatives_;
976 if (constraintDerivativesOld_.empty())
978 constraintDerivativesOld_.setSize(constraintDerivatives_.size());
980 forAll(constraintDerivativesOld_, cI)
982 constraintDerivativesOld_[cI] = constraintDerivatives_[cI];
984 correctionOld_ = correction_;
1005 const label nConstraints,
1012 doAllocateLagrangeMultipliers_(true),
1013 includeBoundConstraints_
1015 designVars->upperBounds() && designVars->lowerBounds()
1019 coeffsDict(
type).getOrDefault<bool>(
"includeExtraVars", false)
1021 p_(activeDesignVars_.size(),
Zero),
1022 gs_(nConstraints_, 1),
1025 includeBoundConstraints_ &&
found(
"lTilda") ?
1026 new
scalarField(
"lTilda", *this, activeDesignVars_.size()) :
1032 includeBoundConstraints_ &&
found(
"uTilda") ?
1033 new
scalarField(
"uTilda", *this, activeDesignVars_.size()) :
1037 extraVars_(nullptr),
1039 c_(coeffsDict(
type).getOrDefault<scalar>(
"c", 100)),
1040 deltaP_(activeDesignVars_.size(),
Zero),
1041 deltaLamda_(nConstraints_,
Zero),
1042 deltaGs_(nConstraints_,
Zero),
1043 deltaLTilda_(nullptr),
1045 deltaUTilda_(nullptr),
1047 deltaExtraVars_(nullptr),
1050 epsMin_(coeffsDict(
type).getOrDefault<scalar>(
"epsMin", 1.
e-07)),
1051 maxNewtonIters_(coeffsDict(
type).getOrDefault<label>(
"maxIters", 1000)),
1054 coeffsDict(
type).getOrDefault<label>(
"maxLineSearchIters", 10)
1056 maxDxIters_(coeffsDict(
type).getOrDefault<label>(
"maxDpIters", 1000)),
1057 relTol_(coeffsDict(
type).getOrDefault<scalar>(
"relTol", 0.01)),
1060 preconditionerNames.getOrDefault
1062 "preconditioner", coeffsDict(
type), preconditioners::
diag 1066 (coeffsDict(
type).getOrDefault<scalar>(
"targetConstraintReduction", 1)),
1067 meritFunctionFile_(nullptr)
1069 Info<<
"Preconditioner type of the SQP subproblem is ::" 1090 LagrangianDerivatives_ = objectiveDerivatives_;
1113 scalar
L = objectiveValue_ + mu_*
sum(
pos(cValues_)*cValues_);
1122 globalSum(objectiveDerivatives_*correction_)
1123 - mu_*
sum(
pos(cValues_)*cValues_);
1132 correctionOld_ = oldCorrection;
1138 if (includeBoundConstraints_)
1140 uTilda_().writeEntry(
"uTilda",
os);
1141 lTilda_().writeEntry(
"lTilda",
os);
tmp< scalarField > preconVectorProduct(const scalarField &vector, autoPtr< scalarField > &precon)
Preconditioner-vector product for CG.
tmp< scalarField > resFuTilda()
Residual of the complementarity slackness for the upper bound constraints.
tmp< scalarField > diagPreconditioner()
Diagonal preconditioner of the deltaP eqn.
void size(const label n)
Older name for setAddressableSize.
virtual void updateOldCorrection(const scalarField &oldCorrection)
Update old correction. Useful for quasi-Newton methods coupled with line search.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
scalarField p_
The set of design variables being updated during the subproblem.
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"))
Type gMin(const FieldField< Field, Type > &f)
void updateYS()
Update the vectors accosiated with the Hessian matrix.
void initialize()
Allocate fields related to constraints.
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
virtual scalar computeMeritFunction()
Compute merit function. Could be different than the objective in the presence of constraints.
autoPtr< scalarField > deltaUTilda_
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual bool addToFile(Ostream &os) const
Write continuation info.
tmp< scalarField > ShermanMorrisonRank1Update(const PtrList< scalarField > &r1Updates, const scalarField &p, const refPtr< scalarField > &diag, const scalarField &mult, label n)
Recursive (naive) implementation of the rank-1 update.
tmp< scalarField > resFlTilda()
Residual of the complementarity slackness for the lower bound constraints.
autoPtr< scalarField > uTilda_
Lagrange multipliers for the upper bound constraints.
bool includeBoundConstraints_
Are bound constraints included?
autoPtr< scalarField > lTilda_
Lagrange multipliers for the lower bound constraints.
#define addProfiling(Name,...)
Define profiling trigger with specified name and description string. The description is generated by ...
virtual bool writeData(Ostream &os) const
Write useful quantities to files.
virtual scalar meritFunctionDirectionalDerivative()
Derivative of the merit function. Could be different than the objective derivative in the presence of...
label preconType_
Which preconditioner to use for the solution of the subproblem.
tmp< scalarField > computeRHSForDeltaX(const scalarField &FDx)
Compute the RHS for the deltaX equation.
constexpr char tab
The tab '\t' character(0x09)
void computeCorrection()
Compute design variables correction.
autoPtr< scalarField > deltaLs_
A class for managing references or pointers (no reference counting)
Macros for easy insertion into run-time selection tables.
void adjustStep(scalar &step, const scalar value, const scalar update)
Adjust step to satisfy cireteria.
virtual scalar meritFunctionConstraintPart() const
Get the part the merit function that depends on the constraints.
void CGforDeltaP(const scalarField &FDx)
CG algorithm for the solution of the deltaP eqn.
#define forAll(list, i)
Loop across all elements in list.
void reset(T *p=nullptr) noexcept
Delete managed object and set to new given pointer.
tmp< scalarField > resFExtraVars()
Residual of the Lagrangian derivative wrt the extra variables.
tmp< scalarField > resFL()
Residual of the gradient of the Lagrangian.
void updateSizes()
Update sizes of fields related to the active design variables.
tmp< scalarField > resFlamda()
Residual of the complementarity slackness for the inequality constraints.
dimensionedScalar pos(const dimensionedScalar &ds)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
const labelList & activeDesignVars_
Map to active design variables.
autoPtr< scalarField > deltaLTilda_
virtual void updateOldCorrection(const scalarField &oldCorrection)
Update old correction. Useful for quasi-newton methods coupled with line search.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
void solveDeltaPEqn()
Solve the equation for deltaX, which is the expensive part of the Newtopn step.
Type gSum(const FieldField< Field, Type > &f)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
void computeNewtonDirection()
Compute direction for the Newton problem.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
PtrList< scalarField > s_
The previous correction. Holds nPrevSteps_ fields.
tmp< scalarField > resFls()
Residual of the lower bounds slackness.
scalar lineSearch()
Perform line search and return max residual corresponding to the updated solution.
tmp< scalarField > computeResiduals()
Compute and return residuals based on the current solution.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
label nPrevSteps_
Number of old corrections and grad differences kept.
bool useSDamping_
Use damping for s to ensure positive-definitiveness.
Base class for Sequantial Quadratic Programming (SQP) methods.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
virtual bool writeData(Ostream &os) const
Write useful quantities to files.
void allocateBoundMultipliers()
Allocate multipliers for the bound constraints.
#define DebugInfo
Report an information message using Foam::Info.
tmp< scalarField > resFz()
Residual of the complementarity slackness for the extra variables.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
const List< word > & names() const noexcept
The list of enum names, in construction order. Same as toc()
Istream and Ostream manipulators taking arguments.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
tmp< scalarField > resFGs()
Residual of the inequality constraint slackness.
OBJstream os(runTime.globalPath()/outputName)
tmp< scalarField > ShermanMorrisonPrecon(const scalarField &vector)
Use the Sherman-Morrison formula to compute the matrix-free product of the approximate inverse of the...
void solveSubproblem()
Solve subproblem using a Newton optimiser.
defineTypeNameAndDebug(combustionModel, 0)
bool useYDamping_
Use damping for s to ensure positive-definitiveness.
An L-BFGS-based SQP algorithm for computing the update of the design variables in the presence of ine...
void updateSolution(const scalar step)
Update the current solution using the known direction and the given step length.
virtual bool writeMeritFunction(const updateMethod &UpdateMethod)
Write info about the merit function.
static const Enum< preconditioners > preconditionerNames
Names of preconditioners for the subproblem.
virtual bool writeAuxiliaryData()
Write merit function information.
autoPtr< scalarField > ls_
Slack variables the lower bound constraints.
#define WarningInFunction
Report a warning using Foam::Warning.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
Mesh data needed to do the Finite Volume discretisation.
PtrList< scalarField > y_
The previous differences of derivatives. Holds nPrevSteps_ fields.
tmp< scalarField > matrixVectorProduct(const scalarField &vector)
Procudt of the LHS of the deltaP eqn with a vector.
The quasi-Newton Limited-memory BFGS formula. Keeps nPrevSteps_ of the y and s vectors and approximat...
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
autoPtr< scalarField > us_
Slack variables the upper bound constraints.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
tmp< scalarField > resFus()
Residual of the upper bounds slackness.
A class for managing temporary objects.
tmp< scalarField > invHFL()
Product of the inverse Hessian with the residual of the gradient of the Lagrangian.
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
void storeOldFields()
Store old fields needed for the next iter.
autoPtr< scalarField > deltaUs_
void allocateLagrangeMultipliers()
Allocate Lagrange multipliers for the inequality constraints.
void zeroUpdates()
Zero the updates computed in the previous optimisation cycle.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
void reset(T *p=nullptr) noexcept
Delete managed pointer and set to new given pointer.
static constexpr const zero Zero
Global zero (0)