47 constrainedOptimisationMethod,
56 void Foam::SQP::allocateMatrices()
82 void Foam::SQP::updateHessian()
87 scalarField LagrangianDerivativesOld = objectiveDerivativesOld_;
88 forAll(constraintDerivatives_, cI)
90 LagrangianDerivatives_ -= lamdas_[cI] * constraintDerivatives_[cI];
91 LagrangianDerivativesOld -= lamdas_[cI] * constraintDerivativesOld_[cI];
93 y.map(LagrangianDerivatives_ - LagrangianDerivativesOld, activeDesignVars_);
94 s.map(correctionOld_, activeDesignVars_);
96 scalar ys = globalSum(
s*
y);
97 if (counter_ == 1 && scaleFirstHessian_)
101 scalar scaleFactor = ys/globalSum(
y*
y);
102 Info<<
"Scaling Hessian with factor " << scaleFactor <<
endl;
103 forAll(activeDesignVars_, varI)
105 HessianOld_[varI][varI] /= scaleFactor;
111 <<
" y*s is negative. Skipping the scaling of the first Hessian"
115 scalar sBs = globalSum(leftMult(
s, HessianOld_)*
s);
119 if (ys < dumpingThreshold_*sBs)
122 <<
" y*s is below threshold. Using damped form" <<
endl;
123 theta = (1 - dumpingThreshold_)*sBs/(sBs - ys);
125 scalarField r(theta*
y + (scalar(1) - theta)*rightMult(HessianOld_,
s));
127 <<
"Unmodified Hessian curvature index " << ys <<
endl;
129 <<
"Modified Hessian curvature index " << globalSum(r*
s) <<
endl;
134 - outerProd(rightMult(HessianOld_,
s), leftMult(
s/sBs, HessianOld_))
135 + outerProd(r, r/globalSum(
s*r));
139 void Foam::SQP::computeLagrangeMultipliersAndCorrect()
141 SquareMatrix<scalar> HessianInv =
inv(Hessian_);
144 Info<<
"Hessian " << Hessian_ <<
endl;
145 Info<<
"HessianInv " << HessianInv <<
endl;
146 label
n = Hessian_.n();
147 SquareMatrix<scalar> test(
n,
Zero);
148 for (label
k = 0;
k <
n;
k++)
150 for (label l = 0; l <
n; l++)
153 for (label i = 0; i <
n; i++)
155 elem += Hessian_[
k][i] * HessianInv[i][l];
160 Info<<
"Validation " << test <<
endl;
164 label nc = constraintDerivatives_.size();
168 activeDerivs.
map(LagrangianDerivatives_, activeDesignVars_);
169 scalarField WgradL = rightMult(HessianInv, activeDerivs);
175 activeConsDerivs.
map(constraintDerivatives_[cI], activeDesignVars_);
176 lamdaRHS[cI] = globalSum(activeConsDerivs * WgradL) - cValues_[cI];
179 Info<<
"lamdaRHS total|deriv part|constraint part "
180 << lamdaRHS[cI] <<
" " << globalSum(activeConsDerivs * WgradL)
181 <<
" " << cValues_[cI] <<
endl;
186 SquareMatrix<scalar> AWA(nc,
Zero);
187 PtrList<scalarField> WA(nc);
188 for (label j = 0; j < nc; j++)
191 gradcJ.
map(constraintDerivatives_[j], activeDesignVars_);
192 WA.set(j,
new scalarField(rightMult(HessianInv, gradcJ)));
193 for (label i = 0; i < nc; i++)
196 gradcI.
map(constraintDerivatives_[i], activeDesignVars_);
197 AWA[i][j] = globalSum(gradcI * WA[j]);
200 SquareMatrix<scalar> invAWA =
inv(AWA);
201 scalarField deltaLamda = rightMult(invAWA, lamdaRHS);
206 Info<<
"lamda update " << deltaLamda <<
endl;
208 lamdas_ += deltaLamda;
214 activeCorrection += WA[cI]*deltaLamda[cI];
216 activeCorrection *= etaHessian_;
219 forAll(activeDesignVars_, varI)
221 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
230 void Foam::SQP::storeOldFields()
232 objectiveDerivativesOld_ = objectiveDerivatives_;
233 if (constraintDerivativesOld_.empty())
235 constraintDerivativesOld_.setSize(constraintDerivatives_.size());
237 forAll(constraintDerivativesOld_, cI)
239 constraintDerivativesOld_[cI] = constraintDerivatives_[cI];
241 correctionOld_ = correction_;
242 HessianOld_ = Hessian_;
246 void Foam::SQP::readFromDict()
248 if (optMethodIODict_.headerOk())
250 optMethodIODict_.readEntry(
"Hessian", Hessian_);
251 optMethodIODict_.readEntry(
"HessianOld", HessianOld_);
252 optMethodIODict_.readEntry
254 "objectiveDerivativesOld",
255 objectiveDerivativesOld_
257 optMethodIODict_.readEntry
259 "constraintDerivativesOld",
260 constraintDerivativesOld_
262 optMethodIODict_.readEntry(
"correctionOld", correctionOld_);
263 optMethodIODict_.readEntry(
"lamdas", lamdas_);
264 optMethodIODict_.readEntry(
"counter", counter_);
265 optMethodIODict_.readEntry(
"eta", eta_);
269 if (activeDesignVars_.empty())
271 activeDesignVars_ =
identity(correction_.size());
285 coeffsDict().getOrDefault<scalar>(
"etaHessian", 1)
287 activeDesignVars_(0),
290 coeffsDict().getOrDefault<
bool>(
"scaleFirstHessian", false)
294 coeffsDict().getOrDefault<scalar>(
"dumpingThreshold", 0.2)
296 LagrangianDerivatives_(0),
299 objectiveDerivativesOld_(0),
300 constraintDerivativesOld_(0),
306 mesh_.time().globalPath()/
"optimisation"/
"objective"/
309 meritFunctionFile_(nullptr),
313 coeffsDict().getOrDefault<scalar>(
"delta", 0.1)
323 Info<<
"\t Did not find explicit definition of active design "
324 <<
"variables. Treating all available ones as active " <<
endl;
349 LagrangianDerivatives_ = objectiveDerivatives_;
356 computeLagrangeMultipliersAndCorrect();
368 if (mu_ <
max(
mag(lamdas_)) + delta_)
370 mu_ =
max(
mag(lamdas_)) + 2*delta_;
373 Info<<
"Updated mu value to " << mu_ <<
endl;
376 scalar
L = objectiveValue_ + mu_*
sum(
mag(cValues_));
385 globalSum(objectiveDerivatives_*correction_)
395 correctionOld_ = oldCorrection;
407 "objectiveDerivativesOld", objectiveDerivativesOld_,
true
410 add<List<scalarField>>
412 "constraintDerivativesOld", constraintDerivativesOld_, true
414 optMethodIODict_.add<
scalarField>(
"correctionOld", correctionOld_,
true);
415 optMethodIODict_.add<
scalarField>(
"lamdas", lamdas_,
true);
416 optMethodIODict_.add<label>(
"counter", counter_,
true);
421 scalar constraintPart =
sum(
mag(cValues_));
422 scalar merit = objectiveValue_ + mu_*constraintPart;
426 unsigned int constraintsSize = lamdas_.size();
427 constraintsSize = constraintsSize*(width + 1) + 2;
430 if (meritFunctionFile_.empty())
432 meritFunctionFile_.reset
438 <<
setw(1) <<
"#" <<
" "
439 <<
setw(width) <<
"merit" <<
" "
440 <<
setw(width) <<
"J" <<
" "
441 <<
setw(constraintsSize) <<
"lamdas" <<
" "
442 <<
setw(constraintsSize) <<
"constraints" <<
" "
443 <<
setw(width) <<
"mu" <<
" "
444 <<
setw(width) <<
"constraintContr" <<
endl;
449 <<
setw(1) << mesh_.time().value() -1 <<
" "
450 <<
setw(width) << merit <<
" "
451 <<
setw(width) << objectiveValue_ <<
" "
457 <<
setw(width) << lamdas_[cI] <<
setw(1) <<
" ";
459 meritFunctionFile_() <<
setw(3) <<
")(";
463 <<
setw(width) << cValues_[cI] <<
setw(1) <<
" ";
465 meritFunctionFile_() <<
setw(2) <<
") ";
466 meritFunctionFile_() <<
setw(width) << mu_ <<
" ";
467 meritFunctionFile_() <<
setw(width) << constraintPart <<
endl;