1 #include <../src/tao/constrained/impls/almm/almm.h> /*I "petsctao.h" I*/ /*I "petscvec.h" I*/ 2 #include <petsctao.h> 3 #include <petsc/private/petscimpl.h> 4 #include <petsc/private/vecimpl.h> 5 6 static PetscErrorCode TaoALMMCombinePrimal_Private(Tao, Vec, Vec, Vec); 7 static PetscErrorCode TaoALMMCombineDual_Private(Tao, Vec, Vec, Vec); 8 static PetscErrorCode TaoALMMSplitPrimal_Private(Tao, Vec, Vec, Vec); 9 static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao); 10 static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao); 11 static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao); 12 13 static PetscErrorCode TaoSolve_ALMM(Tao tao) { 14 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 15 TaoConvergedReason reason; 16 PetscReal updated; 17 18 PetscFunctionBegin; 19 /* reset initial multiplier/slack guess */ 20 if (!tao->recycle) { 21 if (tao->ineq_constrained) { 22 PetscCall(VecZeroEntries(auglag->Ps)); 23 PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->Px, auglag->Ps, auglag->P)); 24 PetscCall(VecZeroEntries(auglag->Yi)); 25 } 26 if (tao->eq_constrained) PetscCall(VecZeroEntries(auglag->Ye)); 27 } 28 29 /* compute initial nonlinear Lagrangian and its derivatives */ 30 PetscCall((*auglag->sub_obj)(tao)); 31 PetscCall(TaoALMMComputeOptimalityNorms_Private(tao)); 32 /* print initial step and check convergence */ 33 PetscCall(PetscInfo(tao, "Solving with %s formulation\n", TaoALMMTypes[auglag->type])); 34 PetscCall(TaoLogConvergenceHistory(tao, auglag->Lval, auglag->gnorm, auglag->cnorm, tao->ksp_its)); 35 PetscCall(TaoMonitor(tao, tao->niter, auglag->fval, auglag->gnorm, auglag->cnorm, 0.0)); 36 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 37 /* set initial penalty factor and inner solver tolerance */ 38 switch (auglag->type) { 39 case TAO_ALMM_CLASSIC: auglag->mu = auglag->mu0; break; 40 case TAO_ALMM_PHR: 41 auglag->cenorm = 0.0; 42 if (tao->eq_constrained) { PetscCall(VecDot(auglag->Ce, auglag->Ce, &auglag->cenorm)); } 43 auglag->cinorm = 0.0; 44 if (tao->ineq_constrained) { 45 PetscCall(VecCopy(auglag->Ci, auglag->Ciwork)); 46 PetscCall(VecScale(auglag->Ciwork, -1.0)); 47 PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork)); 48 PetscCall(VecDot(auglag->Ciwork, auglag->Ciwork, &auglag->cinorm)); 49 } 50 /* determine initial penalty factor based on the balance of constraint violation and objective function value */ 51 auglag->mu = PetscMax(1.e-6, PetscMin(10.0, 2.0 * PetscAbsReal(auglag->fval) / (auglag->cenorm + auglag->cinorm))); 52 break; 53 default: break; 54 } 55 auglag->gtol = auglag->gtol0; 56 PetscCall(PetscInfo(tao, "Initial penalty: %.2g\n", (double)auglag->mu)); 57 58 /* start aug-lag outer loop */ 59 while (tao->reason == TAO_CONTINUE_ITERATING) { 60 ++tao->niter; 61 /* update subsolver tolerance */ 62 PetscCall(PetscInfo(tao, "Subsolver tolerance: ||G|| <= %e\n", (double)auglag->gtol)); 63 PetscCall(TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0)); 64 /* solve the bound-constrained or unconstrained subproblem */ 65 PetscCall(TaoSolve(auglag->subsolver)); 66 PetscCall(TaoGetConvergedReason(auglag->subsolver, &reason)); 67 tao->ksp_its += auglag->subsolver->ksp_its; 68 if (reason != TAO_CONVERGED_GATOL) { PetscCall(PetscInfo(tao, "Subsolver failed to converge, reason: %s\n", TaoConvergedReasons[reason])); } 69 /* evaluate solution and test convergence */ 70 PetscCall((*auglag->sub_obj)(tao)); 71 PetscCall(TaoALMMComputeOptimalityNorms_Private(tao)); 72 /* decide whether to update multipliers or not */ 73 updated = 0.0; 74 if (auglag->cnorm <= auglag->ytol) { 75 PetscCall(PetscInfo(tao, "Multipliers updated: ||C|| <= %e\n", (double)auglag->ytol)); 76 /* constraints are good, update multipliers and convergence tolerances */ 77 if (tao->eq_constrained) { 78 PetscCall(VecAXPY(auglag->Ye, auglag->mu, auglag->Ce)); 79 PetscCall(VecSet(auglag->Cework, auglag->ye_max)); 80 PetscCall(VecPointwiseMin(auglag->Ye, auglag->Cework, auglag->Ye)); 81 PetscCall(VecSet(auglag->Cework, auglag->ye_min)); 82 PetscCall(VecPointwiseMax(auglag->Ye, auglag->Cework, auglag->Ye)); 83 } 84 if (tao->ineq_constrained) { 85 PetscCall(VecAXPY(auglag->Yi, auglag->mu, auglag->Ci)); 86 PetscCall(VecSet(auglag->Ciwork, auglag->yi_max)); 87 PetscCall(VecPointwiseMin(auglag->Yi, auglag->Ciwork, auglag->Yi)); 88 PetscCall(VecSet(auglag->Ciwork, auglag->yi_min)); 89 PetscCall(VecPointwiseMax(auglag->Yi, auglag->Ciwork, auglag->Yi)); 90 } 91 /* tolerances are updated only for non-PHR methods */ 92 if (auglag->type != TAO_ALMM_PHR) { 93 auglag->ytol = PetscMax(tao->catol, auglag->ytol / PetscPowReal(auglag->mu, auglag->mu_pow_good)); 94 auglag->gtol = PetscMax(tao->gatol, auglag->gtol / auglag->mu); 95 } 96 updated = 1.0; 97 } else { 98 /* constraints are bad, update penalty factor */ 99 auglag->mu = PetscMin(auglag->mu_max, auglag->mu_fac * auglag->mu); 100 /* tolerances are reset only for non-PHR methods */ 101 if (auglag->type != TAO_ALMM_PHR) { 102 auglag->ytol = PetscMax(tao->catol, 0.1 / PetscPowReal(auglag->mu, auglag->mu_pow_bad)); 103 auglag->gtol = PetscMax(tao->gatol, 1.0 / auglag->mu); 104 } 105 PetscCall(PetscInfo(tao, "Penalty increased: mu = %.2g\n", (double)auglag->mu)); 106 } 107 PetscCall(TaoLogConvergenceHistory(tao, auglag->fval, auglag->gnorm, auglag->cnorm, tao->ksp_its)); 108 PetscCall(TaoMonitor(tao, tao->niter, auglag->fval, auglag->gnorm, auglag->cnorm, updated)); 109 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 110 } 111 112 PetscFunctionReturn(0); 113 } 114 115 static PetscErrorCode TaoView_ALMM(Tao tao, PetscViewer viewer) { 116 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 117 PetscBool isascii; 118 119 PetscFunctionBegin; 120 PetscCall(PetscViewerASCIIPushTab(viewer)); 121 PetscCall(TaoView(auglag->subsolver, viewer)); 122 PetscCall(PetscViewerASCIIPopTab(viewer)); 123 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 124 if (isascii) { 125 PetscCall(PetscViewerASCIIPushTab(viewer)); 126 PetscCall(PetscViewerASCIIPrintf(viewer, "ALMM Formulation Type: %s\n", TaoALMMTypes[auglag->type])); 127 PetscCall(PetscViewerASCIIPopTab(viewer)); 128 } 129 PetscFunctionReturn(0); 130 } 131 132 static PetscErrorCode TaoSetUp_ALMM(Tao tao) { 133 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 134 VecType vec_type; 135 Vec SL, SU; 136 PetscBool is_cg = PETSC_FALSE, is_lmvm = PETSC_FALSE; 137 138 PetscFunctionBegin; 139 PetscCheck(!tao->ineq_doublesided, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "TAOALMM does not support double-sided inequality constraint definition. Please restructure your inequality constrainst to fit the form c(x) >= 0."); 140 PetscCheck(tao->eq_constrained || tao->ineq_constrained, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "Equality and/or inequality constraints must be defined before solver setup."); 141 PetscCall(TaoComputeVariableBounds(tao)); 142 /* alias base vectors and create extras */ 143 PetscCall(VecGetType(tao->solution, &vec_type)); 144 auglag->Px = tao->solution; 145 if (!tao->gradient) { /* base gradient */ 146 PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 147 } 148 auglag->LgradX = tao->gradient; 149 if (!auglag->Xwork) { /* opt var work vector */ 150 PetscCall(VecDuplicate(tao->solution, &auglag->Xwork)); 151 } 152 if (tao->eq_constrained) { 153 auglag->Ce = tao->constraints_equality; 154 auglag->Ae = tao->jacobian_equality; 155 if (!auglag->Ye) { /* equality multipliers */ 156 PetscCall(VecDuplicate(auglag->Ce, &auglag->Ye)); 157 } 158 if (!auglag->Cework) { PetscCall(VecDuplicate(auglag->Ce, &auglag->Cework)); } 159 } 160 if (tao->ineq_constrained) { 161 auglag->Ci = tao->constraints_inequality; 162 auglag->Ai = tao->jacobian_inequality; 163 if (!auglag->Yi) { /* inequality multipliers */ 164 PetscCall(VecDuplicate(auglag->Ci, &auglag->Yi)); 165 } 166 if (!auglag->Ciwork) { PetscCall(VecDuplicate(auglag->Ci, &auglag->Ciwork)); } 167 if (!auglag->Cizero) { 168 PetscCall(VecDuplicate(auglag->Ci, &auglag->Cizero)); 169 PetscCall(VecZeroEntries(auglag->Cizero)); 170 } 171 if (!auglag->Ps) { /* slack vars */ 172 PetscCall(VecDuplicate(auglag->Ci, &auglag->Ps)); 173 } 174 if (!auglag->LgradS) { /* slack component of Lagrangian gradient */ 175 PetscCall(VecDuplicate(auglag->Ci, &auglag->LgradS)); 176 } 177 /* create vector for combined primal space and the associated communication objects */ 178 if (!auglag->P) { 179 PetscCall(PetscMalloc1(2, &auglag->Parr)); 180 auglag->Parr[0] = auglag->Px; 181 auglag->Parr[1] = auglag->Ps; 182 PetscCall(VecConcatenate(2, auglag->Parr, &auglag->P, &auglag->Pis)); 183 PetscCall(PetscMalloc1(2, &auglag->Pscatter)); 184 PetscCall(VecScatterCreate(auglag->P, auglag->Pis[0], auglag->Px, NULL, &auglag->Pscatter[0])); 185 PetscCall(VecScatterCreate(auglag->P, auglag->Pis[1], auglag->Ps, NULL, &auglag->Pscatter[1])); 186 } 187 if (tao->eq_constrained) { 188 /* create vector for combined dual space and the associated communication objects */ 189 if (!auglag->Y) { 190 PetscCall(PetscMalloc1(2, &auglag->Yarr)); 191 auglag->Yarr[0] = auglag->Ye; 192 auglag->Yarr[1] = auglag->Yi; 193 PetscCall(VecConcatenate(2, auglag->Yarr, &auglag->Y, &auglag->Yis)); 194 PetscCall(PetscMalloc1(2, &auglag->Yscatter)); 195 PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[0], auglag->Ye, NULL, &auglag->Yscatter[0])); 196 PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[1], auglag->Yi, NULL, &auglag->Yscatter[1])); 197 } 198 if (!auglag->C) { PetscCall(VecDuplicate(auglag->Y, &auglag->C)); } 199 } else { 200 if (!auglag->C) { auglag->C = auglag->Ci; } 201 if (!auglag->Y) { auglag->Y = auglag->Yi; } 202 } 203 } else { 204 if (!auglag->P) { auglag->P = auglag->Px; } 205 if (!auglag->G) { auglag->G = auglag->LgradX; } 206 if (!auglag->C) { auglag->C = auglag->Ce; } 207 if (!auglag->Y) { auglag->Y = auglag->Ye; } 208 } 209 /* initialize parameters */ 210 if (auglag->type == TAO_ALMM_PHR) { 211 auglag->mu_fac = 10.0; 212 auglag->yi_min = 0.0; 213 auglag->ytol0 = 0.5; 214 auglag->gtol0 = tao->gatol; 215 if (tao->gatol_changed && tao->catol_changed) { 216 PetscCall(PetscInfo(tao, "TAOALMM with PHR: different gradient and constraint tolerances are not supported, setting catol = gatol\n")); 217 tao->catol = tao->gatol; 218 } 219 } 220 /* set the Lagrangian formulation type for the subsolver */ 221 switch (auglag->type) { 222 case TAO_ALMM_CLASSIC: auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private; break; 223 case TAO_ALMM_PHR: auglag->sub_obj = TaoALMMComputePHRLagAndGradient_Private; break; 224 default: break; 225 } 226 /* set up the subsolver */ 227 PetscCall(TaoSetSolution(auglag->subsolver, auglag->P)); 228 PetscCall(TaoSetObjective(auglag->subsolver, TaoALMMSubsolverObjective_Private, (void *)auglag)); 229 PetscCall(TaoSetObjectiveAndGradient(auglag->subsolver, NULL, TaoALMMSubsolverObjectiveAndGradient_Private, (void *)auglag)); 230 if (tao->bounded) { 231 /* make sure that the subsolver is a bound-constrained method */ 232 PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg)); 233 PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm)); 234 if (is_cg) { 235 PetscCall(TaoSetType(auglag->subsolver, TAOBNCG)); 236 PetscCall(PetscInfo(tao, "TAOCG detected for bound-constrained problem, switching to TAOBNCG instead.")); 237 } 238 if (is_lmvm) { 239 PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS)); 240 PetscCall(PetscInfo(tao, "TAOLMVM detected for bound-constrained problem, switching to TAOBQNLS instead.")); 241 } 242 /* create lower and upper bound clone vectors for subsolver */ 243 if (!auglag->PL) { PetscCall(VecDuplicate(auglag->P, &auglag->PL)); } 244 if (!auglag->PU) { PetscCall(VecDuplicate(auglag->P, &auglag->PU)); } 245 if (tao->ineq_constrained) { 246 /* create lower and upper bounds for slack, set lower to 0 */ 247 PetscCall(VecDuplicate(auglag->Ci, &SL)); 248 PetscCall(VecSet(SL, 0.0)); 249 PetscCall(VecDuplicate(auglag->Ci, &SU)); 250 PetscCall(VecSet(SU, PETSC_INFINITY)); 251 /* combine opt var bounds with slack bounds */ 252 PetscCall(TaoALMMCombinePrimal_Private(tao, tao->XL, SL, auglag->PL)); 253 PetscCall(TaoALMMCombinePrimal_Private(tao, tao->XU, SU, auglag->PU)); 254 /* destroy work vectors */ 255 PetscCall(VecDestroy(&SL)); 256 PetscCall(VecDestroy(&SU)); 257 } else { 258 /* no inequality constraints, just copy bounds into the subsolver */ 259 PetscCall(VecCopy(tao->XL, auglag->PL)); 260 PetscCall(VecCopy(tao->XU, auglag->PU)); 261 } 262 PetscCall(TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU)); 263 } 264 PetscCall(TaoSetUp(auglag->subsolver)); 265 auglag->G = auglag->subsolver->gradient; 266 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode TaoDestroy_ALMM(Tao tao) { 271 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 272 273 PetscFunctionBegin; 274 PetscCall(TaoDestroy(&auglag->subsolver)); 275 if (tao->setupcalled) { 276 PetscCall(VecDestroy(&auglag->Xwork)); /* opt work */ 277 if (tao->eq_constrained) { 278 PetscCall(VecDestroy(&auglag->Ye)); /* equality multipliers */ 279 PetscCall(VecDestroy(&auglag->Cework)); /* equality work vector */ 280 } 281 if (tao->ineq_constrained) { 282 PetscCall(VecDestroy(&auglag->Ps)); /* slack vars */ 283 auglag->Parr[0] = NULL; /* clear pointer to tao->solution, will be destroyed by TaoDestroy() shell */ 284 PetscCall(PetscFree(auglag->Parr)); /* array of primal vectors */ 285 PetscCall(VecDestroy(&auglag->LgradS)); /* slack grad */ 286 PetscCall(VecDestroy(&auglag->Cizero)); /* zero vector for pointwise max */ 287 PetscCall(VecDestroy(&auglag->Yi)); /* inequality multipliers */ 288 PetscCall(VecDestroy(&auglag->Ciwork)); /* inequality work vector */ 289 PetscCall(VecDestroy(&auglag->P)); /* combo primal */ 290 PetscCall(ISDestroy(&auglag->Pis[0])); /* index set for X inside P */ 291 PetscCall(ISDestroy(&auglag->Pis[1])); /* index set for S inside P */ 292 PetscCall(PetscFree(auglag->Pis)); /* array of P index sets */ 293 PetscCall(VecScatterDestroy(&auglag->Pscatter[0])); 294 PetscCall(VecScatterDestroy(&auglag->Pscatter[1])); 295 PetscCall(PetscFree(auglag->Pscatter)); 296 if (tao->eq_constrained) { 297 PetscCall(VecDestroy(&auglag->Y)); /* combo multipliers */ 298 PetscCall(PetscFree(auglag->Yarr)); /* array of dual vectors */ 299 PetscCall(VecDestroy(&auglag->C)); /* combo constraints */ 300 PetscCall(ISDestroy(&auglag->Yis[0])); /* index set for Ye inside Y */ 301 PetscCall(ISDestroy(&auglag->Yis[1])); /* index set for Yi inside Y */ 302 PetscCall(PetscFree(auglag->Yis)); 303 PetscCall(VecScatterDestroy(&auglag->Yscatter[0])); 304 PetscCall(VecScatterDestroy(&auglag->Yscatter[1])); 305 PetscCall(PetscFree(auglag->Yscatter)); 306 } 307 } 308 if (tao->bounded) { 309 PetscCall(VecDestroy(&auglag->PL)); /* lower bounds for subsolver */ 310 PetscCall(VecDestroy(&auglag->PU)); /* upper bounds for subsolver */ 311 } 312 } 313 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", NULL)); 314 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", NULL)); 315 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", NULL)); 316 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", NULL)); 317 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", NULL)); 318 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", NULL)); 319 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", NULL)); 320 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", NULL)); 321 PetscCall(PetscFree(tao->data)); 322 PetscFunctionReturn(0); 323 } 324 325 static PetscErrorCode TaoSetFromOptions_ALMM(Tao tao, PetscOptionItems *PetscOptionsObject) { 326 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 327 PetscInt i; 328 329 PetscFunctionBegin; 330 PetscOptionsHeadBegin(PetscOptionsObject, "Augmented Lagrangian multiplier method solves problems with general constraints by converting them into a sequence of unconstrained problems."); 331 PetscCall(PetscOptionsReal("-tao_almm_mu_init", "initial penalty parameter", "", auglag->mu0, &auglag->mu0, NULL)); 332 PetscCall(PetscOptionsReal("-tao_almm_mu_factor", "increase factor for the penalty parameter", "", auglag->mu_fac, &auglag->mu_fac, NULL)); 333 PetscCall(PetscOptionsReal("-tao_almm_mu_power_good", "exponential for penalty parameter when multiplier update is accepted", "", auglag->mu_pow_good, &auglag->mu_pow_good, NULL)); 334 PetscCall(PetscOptionsReal("-tao_almm_mu_power_bad", "exponential for penalty parameter when multiplier update is rejected", "", auglag->mu_pow_bad, &auglag->mu_pow_bad, NULL)); 335 PetscCall(PetscOptionsReal("-tao_almm_mu_max", "maximum safeguard for penalty parameter updates", "", auglag->mu_max, &auglag->mu_max, NULL)); 336 PetscCall(PetscOptionsReal("-tao_almm_ye_min", "minimum safeguard for equality multiplier updates", "", auglag->ye_min, &auglag->ye_min, NULL)); 337 PetscCall(PetscOptionsReal("-tao_almm_ye_max", "maximum safeguard for equality multipliers updates", "", auglag->ye_max, &auglag->ye_max, NULL)); 338 PetscCall(PetscOptionsReal("-tao_almm_yi_min", "minimum safeguard for inequality multipliers updates", "", auglag->yi_min, &auglag->yi_min, NULL)); 339 PetscCall(PetscOptionsReal("-tao_almm_yi_max", "maximum safeguard for inequality multipliers updates", "", auglag->yi_max, &auglag->yi_max, NULL)); 340 PetscCall(PetscOptionsEnum("-tao_almm_type", "augmented Lagrangian formulation type for the subproblem", "TaoALMMType", TaoALMMTypes, (PetscEnum)auglag->type, (PetscEnum *)&auglag->type, NULL)); 341 PetscOptionsHeadEnd(); 342 PetscCall(TaoSetOptionsPrefix(auglag->subsolver, ((PetscObject)tao)->prefix)); 343 PetscCall(TaoAppendOptionsPrefix(auglag->subsolver, "tao_almm_subsolver_")); 344 PetscCall(TaoSetFromOptions(auglag->subsolver)); 345 for (i = 0; i < tao->numbermonitors; i++) { 346 PetscCall(PetscObjectReference((PetscObject)tao->monitorcontext[i])); 347 PetscCall(TaoSetMonitor(auglag->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i])); 348 if (tao->monitor[i] == TaoMonitorDefault || tao->monitor[i] == TaoDefaultCMonitor || tao->monitor[i] == TaoDefaultGMonitor || tao->monitor[i] == TaoDefaultSMonitor) { auglag->info = PETSC_TRUE; } 349 } 350 PetscFunctionReturn(0); 351 } 352 353 /* -------------------------------------------------------- */ 354 355 /*MC 356 TaoALMM - Augmented Lagrangian multiplier method for solving nonlinear optimization problems with general constraints. 357 358 Options Database Keys: 359 + -tao_almm_mu_init <real> - initial penalty parameter (default: 10.) 360 . -tao_almm_mu_factor <real> - increase factor for the penalty parameter (default: 100.) 361 . -tao_almm_mu_max <real> - maximum safeguard for penalty parameter updates (default: 1.e20) 362 . -tao_almm_mu_power_good <real> - exponential for penalty parameter when multiplier update is accepted (default: 0.9) 363 . -tao_almm_mu_power_bad <real> - exponential for penalty parameter when multiplier update is rejected (default: 0.1) 364 . -tao_almm_ye_min <real> - minimum safeguard for equality multiplier updates (default: -1.e20) 365 . -tao_almm_ye_max <real> - maximum safeguard for equality multiplier updates (default: 1.e20) 366 . -tao_almm_yi_min <real> - minimum safeguard for inequality multiplier updates (default: -1.e20) 367 . -tao_almm_yi_max <real> - maximum safeguard for inequality multiplier updates (default: 1.e20) 368 - -tao_almm_type <classic,phr> - change formulation of the augmented Lagrangian merit function for the subproblem (default: classic) 369 370 Level: beginner 371 372 Notes: 373 This method converts a constrained problem into a sequence of unconstrained problems via the augmented 374 Lagrangian merit function. Bound constraints are pushed down to the subproblem without any modifications. 375 376 Two formulations are offered for the subproblem: canonical Hestenes-Powell augmented Lagrangian with slack 377 variables for inequality constraints, and a slack-less Powell-Hestenes-Rockafellar (PHR) formulation utilizing a 378 pointwise max() penalty on inequality constraints. The canonical augmented Lagrangian formulation typically 379 converges faster for most problems. However, PHR may be desirable for problems featuring a large number 380 of inequality constraints because it avoids inflating the size of the subproblem with slack variables. 381 382 The subproblem is solved using a nested first-order TAO solver. The user can retrieve a pointer to 383 the subsolver via TaoALMMGetSubsolver() or pass command line arguments to it using the 384 "-tao_almm_subsolver_" prefix. Currently, TaoALMM does not support second-order methods for the 385 subproblem. It is also highly recommended that the subsolver chosen by the user utilize a trust-region 386 strategy for globalization (default: TAOBQNKTR) especially if the outer problem features bound constraints. 387 388 .vb 389 while unconverged 390 solve argmin_x L(x) s.t. l <= x <= u 391 if ||c|| <= y_tol 392 if ||c|| <= c_tol && ||Lgrad|| <= g_tol: 393 problem converged, return solution 394 else 395 constraints sufficiently improved 396 update multipliers and tighten tolerances 397 endif 398 else 399 constraints did not improve 400 update penalty and loosen tolerances 401 endif 402 endwhile 403 .ve 404 405 .seealso: `TaoALMMGetType()`, `TaoALMMSetType()`, `TaoALMMSetSubsolver()`, `TaoALMMGetSubsolver()`, 406 `TaoALMMGetMultipliers()`, `TaoALMMSetMultipliers()`, `TaoALMMGetPrimalIS()`, `TaoALMMGetDualIS()` 407 M*/ 408 PETSC_EXTERN PetscErrorCode TaoCreate_ALMM(Tao tao) { 409 TAO_ALMM *auglag; 410 411 PetscFunctionBegin; 412 PetscCall(PetscNewLog(tao, &auglag)); 413 414 tao->ops->destroy = TaoDestroy_ALMM; 415 tao->ops->setup = TaoSetUp_ALMM; 416 tao->ops->setfromoptions = TaoSetFromOptions_ALMM; 417 tao->ops->view = TaoView_ALMM; 418 tao->ops->solve = TaoSolve_ALMM; 419 420 tao->gatol = 1.e-5; 421 tao->grtol = 0.0; 422 tao->gttol = 0.0; 423 tao->catol = 1.e-5; 424 tao->crtol = 0.0; 425 426 tao->data = (void *)auglag; 427 auglag->parent = tao; 428 auglag->mu0 = 10.0; 429 auglag->mu = auglag->mu0; 430 auglag->mu_fac = 10.0; 431 auglag->mu_max = PETSC_INFINITY; 432 auglag->mu_pow_good = 0.9; 433 auglag->mu_pow_bad = 0.1; 434 auglag->ye_min = PETSC_NINFINITY; 435 auglag->ye_max = PETSC_INFINITY; 436 auglag->yi_min = PETSC_NINFINITY; 437 auglag->yi_max = PETSC_INFINITY; 438 auglag->ytol0 = 0.1 / PetscPowReal(auglag->mu0, auglag->mu_pow_bad); 439 auglag->ytol = auglag->ytol0; 440 auglag->gtol0 = 1.0 / auglag->mu0; 441 auglag->gtol = auglag->gtol0; 442 443 auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private; 444 auglag->type = TAO_ALMM_CLASSIC; 445 auglag->info = PETSC_FALSE; 446 447 PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver)); 448 PetscCall(TaoSetType(auglag->subsolver, TAOBQNKTR)); 449 PetscCall(TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0)); 450 PetscCall(TaoSetMaximumIterations(auglag->subsolver, 1000)); 451 PetscCall(TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000)); 452 PetscCall(TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY)); 453 PetscCall(PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1)); 454 455 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private)); 456 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private)); 457 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private)); 458 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private)); 459 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private)); 460 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private)); 461 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private)); 462 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private)); 463 PetscFunctionReturn(0); 464 } 465 466 static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P) { 467 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 468 469 PetscFunctionBegin; 470 if (tao->ineq_constrained) { 471 PetscCall(VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE)); 472 PetscCall(VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE)); 473 PetscCall(VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE)); 474 PetscCall(VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE)); 475 } else { 476 PetscCall(VecCopy(X, P)); 477 } 478 PetscFunctionReturn(0); 479 } 480 481 static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y) { 482 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 483 484 PetscFunctionBegin; 485 if (tao->eq_constrained) { 486 if (tao->ineq_constrained) { 487 PetscCall(VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE)); 488 PetscCall(VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE)); 489 PetscCall(VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE)); 490 PetscCall(VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE)); 491 } else { 492 PetscCall(VecCopy(EQ, Y)); 493 } 494 } else { 495 PetscCall(VecCopy(IN, Y)); 496 } 497 PetscFunctionReturn(0); 498 } 499 500 static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S) { 501 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 502 503 PetscFunctionBegin; 504 if (tao->ineq_constrained) { 505 PetscCall(VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD)); 506 PetscCall(VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD)); 507 PetscCall(VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD)); 508 PetscCall(VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD)); 509 } else { 510 PetscCall(VecCopy(P, X)); 511 } 512 PetscFunctionReturn(0); 513 } 514 515 /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */ 516 static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao) { 517 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 518 519 PetscFunctionBegin; 520 /* if bounded, project the gradient */ 521 if (tao->bounded) PetscCall(VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX)); 522 if (auglag->type == TAO_ALMM_PHR) { 523 PetscCall(VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm)); 524 auglag->cenorm = 0.0; 525 if (tao->eq_constrained) { PetscCall(VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm)); } 526 auglag->cinorm = 0.0; 527 if (tao->ineq_constrained) { 528 PetscCall(VecCopy(auglag->Yi, auglag->Ciwork)); 529 PetscCall(VecScale(auglag->Ciwork, -1.0 / auglag->mu)); 530 PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork)); 531 PetscCall(VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm)); 532 } 533 auglag->cnorm_old = auglag->cnorm; 534 auglag->cnorm = PetscMax(auglag->cenorm, auglag->cinorm); 535 auglag->ytol = auglag->ytol0 * auglag->cnorm_old; 536 } else { 537 PetscCall(VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm)); 538 PetscCall(VecNorm(auglag->C, NORM_2, &auglag->cnorm)); 539 } 540 PetscFunctionReturn(0); 541 } 542 543 static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao, Vec P) { 544 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 545 546 PetscFunctionBegin; 547 /* split solution into primal and slack components */ 548 PetscCall(TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps)); 549 550 /* compute f, df/dx and the constraints */ 551 PetscCall(TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX)); 552 if (tao->eq_constrained) { 553 PetscCall(TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce)); 554 PetscCall(TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae)); 555 } 556 if (tao->ineq_constrained) { 557 PetscCall(TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci)); 558 PetscCall(TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai)); 559 switch (auglag->type) { 560 case TAO_ALMM_CLASSIC: 561 /* classic formulation converts inequality to equality constraints via slack variables */ 562 PetscCall(VecAXPY(auglag->Ci, -1.0, auglag->Ps)); 563 break; 564 case TAO_ALMM_PHR: 565 /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */ 566 PetscCall(VecScale(auglag->Ci, -1.0)); 567 PetscCall(MatScale(auglag->Ai, -1.0)); 568 break; 569 default: break; 570 } 571 } 572 /* combine constraints into one vector */ 573 PetscCall(TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C)); 574 PetscFunctionReturn(0); 575 } 576 577 /* 578 Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmin(0, Ci + Yi/mu)^T pmin(0, Ci + Yi/mu)] 579 580 dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmin(0, Ci + Yi/mu)^T Ai] 581 582 dLphr/dS = 0 583 */ 584 static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao) { 585 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 586 PetscReal eq_norm = 0.0, ineq_norm = 0.0; 587 588 PetscFunctionBegin; 589 PetscCall(TaoALMMEvaluateIterate_Private(tao, auglag->P)); 590 if (tao->eq_constrained) { 591 /* Ce_work = mu*(Ce + Ye/mu) */ 592 PetscCall(VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce)); 593 PetscCall(VecDot(auglag->Cework, auglag->Cework, &eq_norm)); /* contribution to scalar Lagrangian */ 594 PetscCall(VecScale(auglag->Cework, auglag->mu)); 595 /* dL/dX += mu*(Ce + Ye/mu)^T Ae */ 596 PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX)); 597 } 598 if (tao->ineq_constrained) { 599 /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */ 600 PetscCall(VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci)); 601 PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork)); 602 PetscCall(VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm)); /* contribution to scalar Lagrangian */ 603 /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */ 604 PetscCall(VecScale(auglag->Ciwork, auglag->mu)); 605 PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX)); 606 /* dL/dS = 0 because there are no slacks in PHR */ 607 PetscCall(VecZeroEntries(auglag->LgradS)); 608 } 609 /* combine gradient together */ 610 PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G)); 611 /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */ 612 auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm); 613 PetscFunctionReturn(0); 614 } 615 616 /* 617 Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)] 618 619 dLc/dX = dF/dX + Ye^TAe + Yi^TAi + 0.5*mu*[Ce^TAe + (Ci - S)^TAi] 620 621 dLc/dS = -[Yi + mu*(Ci - S)] 622 */ 623 static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao) { 624 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 625 PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0; 626 627 PetscFunctionBegin; 628 PetscCall(TaoALMMEvaluateIterate_Private(tao, auglag->P)); 629 if (tao->eq_constrained) { 630 /* compute scalar contributions */ 631 PetscCall(VecDot(auglag->Ye, auglag->Ce, &yeTce)); 632 PetscCall(VecDot(auglag->Ce, auglag->Ce, &ceTce)); 633 /* dL/dX += ye^T Ae */ 634 PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX)); 635 /* dL/dX += mu * ce^T Ae */ 636 PetscCall(MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork)); 637 PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork)); 638 } 639 if (tao->ineq_constrained) { 640 /* compute scalar contributions */ 641 PetscCall(VecDot(auglag->Yi, auglag->Ci, &yiTcims)); 642 PetscCall(VecDot(auglag->Ci, auglag->Ci, &cimsTcims)); 643 /* dL/dX += yi^T Ai */ 644 PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX)); 645 /* dL/dX += mu * (ci - s)^T Ai */ 646 PetscCall(MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork)); 647 PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork)); 648 /* dL/dS = -[yi + mu*(ci - s)] */ 649 PetscCall(VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi)); 650 PetscCall(VecScale(auglag->LgradS, -1.0)); 651 } 652 /* combine gradient together */ 653 PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G)); 654 /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */ 655 auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims); 656 PetscFunctionReturn(0); 657 } 658 659 PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, void *ctx) { 660 TAO_ALMM *auglag = (TAO_ALMM *)ctx; 661 662 PetscFunctionBegin; 663 PetscCall(VecCopy(P, auglag->P)); 664 PetscCall((*auglag->sub_obj)(auglag->parent)); 665 *Lval = auglag->Lval; 666 PetscFunctionReturn(0); 667 } 668 669 PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx) { 670 TAO_ALMM *auglag = (TAO_ALMM *)ctx; 671 672 PetscFunctionBegin; 673 PetscCall(VecCopy(P, auglag->P)); 674 PetscCall((*auglag->sub_obj)(auglag->parent)); 675 PetscCall(VecCopy(auglag->G, G)); 676 *Lval = auglag->Lval; 677 PetscFunctionReturn(0); 678 } 679