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