47 constrainedOptimisationMethod,
56 void Foam::SQP::allocateMatrices()
86 void Foam::SQP::updateHessian()
91 scalarField LagrangianDerivativesOld = objectiveDerivativesOld_;
92 forAll(constraintDerivatives_, cI)
94 LagrangianDerivatives_ -= lamdas_[cI] * constraintDerivatives_[cI];
95 LagrangianDerivativesOld -= lamdas_[cI] * constraintDerivativesOld_[cI];
97 y.map(LagrangianDerivatives_ - LagrangianDerivativesOld, activeDesignVars_);
98 s.map(correctionOld_, activeDesignVars_);
100 scalar ys = globalSum(
s*
y);
101 if (counter_ == 1 && scaleFirstHessian_)
105 scalar scaleFactor = ys/globalSum(
y*
y);
106 Info<<
"Scaling Hessian with factor " << scaleFactor <<
endl;
107 forAll(activeDesignVars_, varI)
109 HessianOld_[varI][varI] /= scaleFactor;
115 <<
" y*s is negative. Skipping the scaling of the first Hessian"
119 scalar sBs = globalSum(leftMult(
s, HessianOld_)*
s);
123 if (ys < dumpingThreshold_*sBs)
126 <<
" y*s is below threshold. Using damped form" <<
endl;
127 theta = (1 - dumpingThreshold_)*sBs/(sBs - ys);
129 scalarField r(theta*
y + (scalar(1) - theta)*rightMult(HessianOld_,
s));
131 <<
"Unmodified Hessian curvature index " << ys <<
endl;
133 <<
"Modified Hessian curvature index " << globalSum(r*
s) <<
endl;
138 - outerProd(rightMult(HessianOld_,
s), leftMult(
s/sBs, HessianOld_))
139 + outerProd(r, r/globalSum(
s*r));
143 void Foam::SQP::computeLagrangeMultipliersAndCorrect()
145 SquareMatrix<scalar> HessianInv =
inv(Hessian_);
148 Info<<
"Hessian " << Hessian_ <<
endl;
149 Info<<
"HessianInv " << HessianInv <<
endl;
151 SquareMatrix<scalar> test(
n,
Zero);
154 for (
label l = 0; l <
n; l++)
157 for (
label i = 0; i <
n; i++)
159 elem += Hessian_[
k][i] * HessianInv[i][l];
164 Info<<
"Validation " << test <<
endl;
168 label nc = constraintDerivatives_.size();
172 activeDerivs.
map(LagrangianDerivatives_, activeDesignVars_);
173 scalarField WgradL = rightMult(HessianInv, activeDerivs);
179 activeConsDerivs.
map(constraintDerivatives_[cI], activeDesignVars_);
180 lamdaRHS[cI] = globalSum(activeConsDerivs * WgradL) - cValues_[cI];
183 Info<<
"lamdaRHS total|deriv part|constraint part "
184 << lamdaRHS[cI] <<
" " << globalSum(activeConsDerivs * WgradL)
185 <<
" " << cValues_[cI] <<
endl;
190 SquareMatrix<scalar> AWA(nc,
Zero);
191 PtrList<scalarField> WA(nc);
192 for (
label j = 0; j < nc; j++)
195 gradcJ.
map(constraintDerivatives_[j], activeDesignVars_);
196 WA.set(j,
new scalarField(rightMult(HessianInv, gradcJ)));
197 for (
label i = 0; i < nc; i++)
200 gradcI.
map(constraintDerivatives_[i], activeDesignVars_);
201 AWA[i][j] = globalSum(gradcI * WA[j]);
204 SquareMatrix<scalar> invAWA =
inv(AWA);
205 scalarField deltaLamda = rightMult(invAWA, lamdaRHS);
210 Info<<
"lamda update " << deltaLamda <<
endl;
212 lamdas_ += deltaLamda;
218 activeCorrection += WA[cI]*deltaLamda[cI];
220 activeCorrection *= etaHessian_;
223 forAll(activeDesignVars_, varI)
225 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
234 void Foam::SQP::storeOldFields()
236 objectiveDerivativesOld_ = objectiveDerivatives_;
237 if (constraintDerivativesOld_.empty())
239 constraintDerivativesOld_.setSize(constraintDerivatives_.size());
241 forAll(constraintDerivativesOld_, cI)
243 constraintDerivativesOld_[cI] = constraintDerivatives_[cI];
245 correctionOld_ = correction_;
246 HessianOld_ = Hessian_;
250 void Foam::SQP::readFromDict()
252 if (optMethodIODict_.headerOk())
254 optMethodIODict_.readEntry(
"Hessian", Hessian_);
255 optMethodIODict_.readEntry(
"HessianOld", HessianOld_);
256 optMethodIODict_.readEntry
258 "objectiveDerivativesOld",
259 objectiveDerivativesOld_
261 optMethodIODict_.readEntry
263 "constraintDerivativesOld",
264 constraintDerivativesOld_
266 optMethodIODict_.readEntry(
"correctionOld", correctionOld_);
267 optMethodIODict_.readEntry(
"lamdas", lamdas_);
268 optMethodIODict_.readEntry(
"counter", counter_);
269 optMethodIODict_.readEntry(
"eta", eta_);
284 coeffsDict().lookupOrDefault<scalar>(
"etaHessian", 1)
286 activeDesignVars_(0),
289 coeffsDict().lookupOrDefault<
bool>(
"scaleFirstHessian", false)
293 coeffsDict().lookupOrDefault<scalar>(
"dumpingThreshold", 0.2)
295 LagrangianDerivatives_(0),
298 objectiveDerivativesOld_(0),
299 constraintDerivativesOld_(0),
305 mesh_.time().globalPath()/
"optimisation"/
"objective"/
308 meritFunctionFile_(nullptr),
312 coeffsDict().lookupOrDefault<scalar>(
"delta", 0.1)
322 Info<<
"\t Did not find explicit definition of active design "
323 <<
"variables. Treating all available ones as active " <<
endl;
348 LagrangianDerivatives_ = objectiveDerivatives_;
355 computeLagrangeMultipliersAndCorrect();
367 if (mu_ <
max(
mag(lamdas_)) + delta_)
369 mu_ =
max(
mag(lamdas_)) + 2*delta_;
372 Info<<
"Updated mu value to " << mu_ <<
endl;
375 scalar
L = objectiveValue_ + mu_*
sum(
mag(cValues_));
384 globalSum(objectiveDerivatives_*correction_)
394 correctionOld_ = oldCorrection;
406 "objectiveDerivativesOld", objectiveDerivativesOld_,
true
409 add<List<scalarField>>
411 "constraintDerivativesOld", constraintDerivativesOld_, true
413 optMethodIODict_.add<
scalarField>(
"correctionOld", correctionOld_,
true);
414 optMethodIODict_.add<
scalarField>(
"lamdas", lamdas_,
true);
415 optMethodIODict_.add<
label>(
"counter", counter_,
true);
420 scalar constraintPart =
sum(
mag(cValues_));
421 scalar merit = objectiveValue_ + mu_*constraintPart;
425 unsigned int constraintsSize = lamdas_.size();
426 constraintsSize = constraintsSize*(width + 1) + 2;
429 if (meritFunctionFile_.empty())
431 meritFunctionFile_.reset
437 <<
setw(1) <<
"#" <<
" "
438 <<
setw(width) <<
"merit" <<
" "
439 <<
setw(width) <<
"J" <<
" "
440 <<
setw(constraintsSize) <<
"lamdas" <<
" "
441 <<
setw(constraintsSize) <<
"constraints" <<
" "
442 <<
setw(width) <<
"mu" <<
" "
443 <<
setw(width) <<
"constraintContr" <<
endl;
448 <<
setw(1) << mesh_.time().value() -1 <<
" "
449 <<
setw(width) << merit <<
" "
450 <<
setw(width) << objectiveValue_ <<
" "
456 <<
setw(width) << lamdas_[cI] <<
setw(1) <<
" ";
458 meritFunctionFile_() <<
setw(3) <<
")(";
462 <<
setw(width) << cValues_[cI] <<
setw(1) <<
" ";
464 meritFunctionFile_() <<
setw(2) <<
") ";
465 meritFunctionFile_() <<
setw(width) << mu_ <<
" ";
466 meritFunctionFile_() <<
setw(width) << constraintPart <<
endl;