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