Actual source code: almm.c

  1: #include <../src/tao/constrained/impls/almm/almm.h>
  2: #include <petsctao.h>
  3: #include <petsc/private/petscimpl.h>
  4: #include <petsc/private/vecimpl.h>

  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);

 13: static PetscErrorCode TaoSolve_ALMM(Tao tao)
 14: {
 15:   TAO_ALMM          *auglag = (TAO_ALMM *)tao->data;
 16:   TaoConvergedReason reason;
 17:   PetscReal          updated;

 19:   /* reset initial multiplier/slack guess */
 20:   if (!tao->recycle) {
 21:     if (tao->ineq_constrained) {
 22:       VecZeroEntries(auglag->Ps);
 23:       TaoALMMCombinePrimal_Private(tao, auglag->Px, auglag->Ps, auglag->P);
 24:       VecSet(auglag->Yi, 1.0);
 25:     }
 26:     if (tao->eq_constrained) VecSet(auglag->Ye, 1.0);
 27:   }

 29:   /* compute initial nonlinear Lagrangian and its derivatives */
 30:   (*auglag->sub_obj)(tao);
 31:   TaoALMMComputeOptimalityNorms_Private(tao);
 32:   /* print initial step and check convergence */
 33:   PetscInfo(tao, "Solving with %s formulation\n", TaoALMMTypes[auglag->type]);
 34:   TaoLogConvergenceHistory(tao, auglag->Lval, auglag->gnorm, auglag->cnorm, tao->ksp_its);
 35:   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:
 40:     auglag->mu = auglag->mu0;
 41:     break;
 42:   case TAO_ALMM_PHR:
 43:     auglag->cenorm = 0.0;
 44:     if (tao->eq_constrained) VecDot(auglag->Ce, auglag->Ce, &auglag->cenorm);
 45:     auglag->cinorm = 0.0;
 46:     if (tao->ineq_constrained) {
 47:       VecCopy(auglag->Ci, auglag->Ciwork);
 48:       VecScale(auglag->Ciwork, -1.0);
 49:       VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork);
 50:       VecDot(auglag->Ciwork, auglag->Ciwork, &auglag->cinorm);
 51:     }
 52:     /* determine initial penalty factor based on the balance of constraint violation and objective function value */
 53:     auglag->mu = PetscMax(1.e-6, PetscMin(10.0, 2.0 * PetscAbsReal(auglag->fval) / (auglag->cenorm + auglag->cinorm)));
 54:     break;
 55:   default:
 56:     break;
 57:   }
 58:   auglag->gtol = auglag->gtol0;
 59:   PetscInfo(tao, "Initial penalty: %.2g\n", (double)auglag->mu);

 61:   /* start aug-lag outer loop */
 62:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 63:     ++tao->niter;
 64:     /* update subsolver tolerance */
 65:     PetscInfo(tao, "Subsolver tolerance: ||G|| <= %e\n", (double)auglag->gtol);
 66:     TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0);
 67:     /* solve the bound-constrained or unconstrained subproblem */
 68:     VecCopy(auglag->P, auglag->subsolver->solution);
 69:     TaoSolve(auglag->subsolver);
 70:     VecCopy(auglag->subsolver->solution, auglag->P);
 71:     TaoGetConvergedReason(auglag->subsolver, &reason);
 72:     tao->ksp_its += auglag->subsolver->ksp_its;
 73:     if (reason != TAO_CONVERGED_GATOL) PetscInfo(tao, "Subsolver failed to converge, reason: %s\n", TaoConvergedReasons[reason]);
 74:     /* evaluate solution and test convergence */
 75:     (*auglag->sub_obj)(tao);
 76:     TaoALMMComputeOptimalityNorms_Private(tao);
 77:     /* decide whether to update multipliers or not */
 78:     updated = 0.0;
 79:     if (auglag->cnorm <= auglag->ytol) {
 80:       PetscInfo(tao, "Multipliers updated: ||C|| <= %e\n", (double)auglag->ytol);
 81:       /* constraints are good, update multipliers and convergence tolerances */
 82:       if (tao->eq_constrained) {
 83:         VecAXPY(auglag->Ye, auglag->mu, auglag->Ce);
 84:         VecSet(auglag->Cework, auglag->ye_max);
 85:         VecPointwiseMin(auglag->Ye, auglag->Cework, auglag->Ye);
 86:         VecSet(auglag->Cework, auglag->ye_min);
 87:         VecPointwiseMax(auglag->Ye, auglag->Cework, auglag->Ye);
 88:       }
 89:       if (tao->ineq_constrained) {
 90:         VecAXPY(auglag->Yi, auglag->mu, auglag->Ci);
 91:         VecSet(auglag->Ciwork, auglag->yi_max);
 92:         VecPointwiseMin(auglag->Yi, auglag->Ciwork, auglag->Yi);
 93:         VecSet(auglag->Ciwork, auglag->yi_min);
 94:         VecPointwiseMax(auglag->Yi, auglag->Ciwork, auglag->Yi);
 95:       }
 96:       /* tolerances are updated only for non-PHR methods */
 97:       if (auglag->type != TAO_ALMM_PHR) {
 98:         auglag->ytol = PetscMax(tao->catol, auglag->ytol / PetscPowReal(auglag->mu, auglag->mu_pow_good));
 99:         auglag->gtol = PetscMax(tao->gatol, auglag->gtol / auglag->mu);
100:       }
101:       updated = 1.0;
102:     } else {
103:       /* constraints are bad, update penalty factor */
104:       auglag->mu = PetscMin(auglag->mu_max, auglag->mu_fac * auglag->mu);
105:       /* tolerances are reset only for non-PHR methods */
106:       if (auglag->type != TAO_ALMM_PHR) {
107:         auglag->ytol = PetscMax(tao->catol, 0.1 / PetscPowReal(auglag->mu, auglag->mu_pow_bad));
108:         auglag->gtol = PetscMax(tao->gatol, 1.0 / auglag->mu);
109:       }
110:       PetscInfo(tao, "Penalty increased: mu = %.2g\n", (double)auglag->mu);
111:     }
112:     TaoLogConvergenceHistory(tao, auglag->fval, auglag->gnorm, auglag->cnorm, tao->ksp_its);
113:     TaoMonitor(tao, tao->niter, auglag->fval, auglag->gnorm, auglag->cnorm, updated);
114:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
115:   }

117:   return 0;
118: }

120: static PetscErrorCode TaoView_ALMM(Tao tao, PetscViewer viewer)
121: {
122:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
123:   PetscBool isascii;

125:   PetscViewerASCIIPushTab(viewer);
126:   TaoView(auglag->subsolver, viewer);
127:   PetscViewerASCIIPopTab(viewer);
128:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
129:   if (isascii) {
130:     PetscViewerASCIIPushTab(viewer);
131:     PetscViewerASCIIPrintf(viewer, "ALMM Formulation Type: %s\n", TaoALMMTypes[auglag->type]);
132:     PetscViewerASCIIPopTab(viewer);
133:   }
134:   return 0;
135: }

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;

146:   TaoComputeVariableBounds(tao);
147:   /* alias base vectors and create extras */
148:   VecGetType(tao->solution, &vec_type);
149:   auglag->Px = tao->solution;
150:   if (!tao->gradient) { /* base gradient */
151:     VecDuplicate(tao->solution, &tao->gradient);
152:   }
153:   auglag->LgradX = tao->gradient;
154:   if (!auglag->Xwork) { /* opt var work vector */
155:     VecDuplicate(tao->solution, &auglag->Xwork);
156:   }
157:   if (tao->eq_constrained) {
158:     auglag->Ce = tao->constraints_equality;
159:     auglag->Ae = tao->jacobian_equality;
160:     if (!auglag->Ye) { /* equality multipliers */
161:       VecDuplicate(auglag->Ce, &auglag->Ye);
162:     }
163:     if (!auglag->Cework) VecDuplicate(auglag->Ce, &auglag->Cework);
164:   }
165:   if (tao->ineq_constrained) {
166:     auglag->Ci = tao->constraints_inequality;
167:     auglag->Ai = tao->jacobian_inequality;
168:     if (!auglag->Yi) { /* inequality multipliers */
169:       VecDuplicate(auglag->Ci, &auglag->Yi);
170:     }
171:     if (!auglag->Ciwork) VecDuplicate(auglag->Ci, &auglag->Ciwork);
172:     if (!auglag->Cizero) {
173:       VecDuplicate(auglag->Ci, &auglag->Cizero);
174:       VecZeroEntries(auglag->Cizero);
175:     }
176:     if (!auglag->Ps) { /* slack vars */
177:       VecDuplicate(auglag->Ci, &auglag->Ps);
178:     }
179:     if (!auglag->LgradS) { /* slack component of Lagrangian gradient */
180:       VecDuplicate(auglag->Ci, &auglag->LgradS);
181:     }
182:     /* create vector for combined primal space and the associated communication objects */
183:     if (!auglag->P) {
184:       PetscMalloc1(2, &auglag->Parr);
185:       auglag->Parr[0] = auglag->Px;
186:       auglag->Parr[1] = auglag->Ps;
187:       VecConcatenate(2, auglag->Parr, &auglag->P, &auglag->Pis);
188:       PetscMalloc1(2, &auglag->Pscatter);
189:       VecScatterCreate(auglag->P, auglag->Pis[0], auglag->Px, NULL, &auglag->Pscatter[0]);
190:       VecScatterCreate(auglag->P, auglag->Pis[1], auglag->Ps, NULL, &auglag->Pscatter[1]);
191:     }
192:     if (tao->eq_constrained) {
193:       /* create vector for combined dual space and the associated communication objects */
194:       if (!auglag->Y) {
195:         PetscMalloc1(2, &auglag->Yarr);
196:         auglag->Yarr[0] = auglag->Ye;
197:         auglag->Yarr[1] = auglag->Yi;
198:         VecConcatenate(2, auglag->Yarr, &auglag->Y, &auglag->Yis);
199:         PetscMalloc1(2, &auglag->Yscatter);
200:         VecScatterCreate(auglag->Y, auglag->Yis[0], auglag->Ye, NULL, &auglag->Yscatter[0]);
201:         VecScatterCreate(auglag->Y, auglag->Yis[1], auglag->Yi, NULL, &auglag->Yscatter[1]);
202:       }
203:       if (!auglag->C) VecDuplicate(auglag->Y, &auglag->C);
204:     } else {
205:       if (!auglag->C) auglag->C = auglag->Ci;
206:       if (!auglag->Y) auglag->Y = auglag->Yi;
207:     }
208:   } else {
209:     if (!auglag->P) auglag->P = auglag->Px;
210:     if (!auglag->G) auglag->G = auglag->LgradX;
211:     if (!auglag->C) auglag->C = auglag->Ce;
212:     if (!auglag->Y) auglag->Y = auglag->Ye;
213:   }
214:   /* create tao primal solution and gradient to interface with subsolver */
215:   VecDuplicate(auglag->P, &auglag->Psub);
216:   VecDuplicate(auglag->P, &auglag->G);
217:   /* initialize parameters */
218:   if (auglag->type == TAO_ALMM_PHR) {
219:     auglag->mu_fac = 10.0;
220:     auglag->yi_min = 0.0;
221:     auglag->ytol0  = 0.5;
222:     auglag->gtol0  = tao->gatol;
223:     if (tao->gatol_changed && tao->catol_changed) {
224:       PetscInfo(tao, "TAOALMM with PHR: different gradient and constraint tolerances are not supported, setting catol = gatol\n");
225:       tao->catol = tao->gatol;
226:     }
227:   }
228:   /* set the Lagrangian formulation type for the subsolver */
229:   switch (auglag->type) {
230:   case TAO_ALMM_CLASSIC:
231:     auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
232:     break;
233:   case TAO_ALMM_PHR:
234:     auglag->sub_obj = TaoALMMComputePHRLagAndGradient_Private;
235:     break;
236:   default:
237:     break;
238:   }
239:   /* set up the subsolver */
240:   TaoSetSolution(auglag->subsolver, auglag->Psub);
241:   TaoSetObjective(auglag->subsolver, TaoALMMSubsolverObjective_Private, (void *)auglag);
242:   TaoSetObjectiveAndGradient(auglag->subsolver, NULL, TaoALMMSubsolverObjectiveAndGradient_Private, (void *)auglag);
243:   if (tao->bounded) {
244:     /* make sure that the subsolver is a bound-constrained method */
245:     PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg);
246:     PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm);
247:     if (is_cg) {
248:       TaoSetType(auglag->subsolver, TAOBNCG);
249:       PetscInfo(tao, "TAOCG detected for bound-constrained problem, switching to TAOBNCG instead.");
250:     }
251:     if (is_lmvm) {
252:       TaoSetType(auglag->subsolver, TAOBQNLS);
253:       PetscInfo(tao, "TAOLMVM detected for bound-constrained problem, switching to TAOBQNLS instead.");
254:     }
255:     /* create lower and upper bound clone vectors for subsolver */
256:     if (!auglag->PL) VecDuplicate(auglag->P, &auglag->PL);
257:     if (!auglag->PU) VecDuplicate(auglag->P, &auglag->PU);
258:     if (tao->ineq_constrained) {
259:       /* create lower and upper bounds for slack, set lower to 0 */
260:       VecDuplicate(auglag->Ci, &SL);
261:       VecSet(SL, 0.0);
262:       VecDuplicate(auglag->Ci, &SU);
263:       VecSet(SU, PETSC_INFINITY);
264:       /* combine opt var bounds with slack bounds */
265:       TaoALMMCombinePrimal_Private(tao, tao->XL, SL, auglag->PL);
266:       TaoALMMCombinePrimal_Private(tao, tao->XU, SU, auglag->PU);
267:       /* destroy work vectors */
268:       VecDestroy(&SL);
269:       VecDestroy(&SU);
270:     } else {
271:       /* no inequality constraints, just copy bounds into the subsolver */
272:       VecCopy(tao->XL, auglag->PL);
273:       VecCopy(tao->XU, auglag->PU);
274:     }
275:     TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU);
276:   }
277:   TaoSetUp(auglag->subsolver);

279:   return 0;
280: }

282: static PetscErrorCode TaoDestroy_ALMM(Tao tao)
283: {
284:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

286:   TaoDestroy(&auglag->subsolver);
287:   VecDestroy(&auglag->Psub);
288:   VecDestroy(&auglag->G);
289:   if (tao->setupcalled) {
290:     VecDestroy(&auglag->Xwork);
291:     if (tao->eq_constrained) {
292:       VecDestroy(&auglag->Ye);     /* equality multipliers */
293:       VecDestroy(&auglag->Cework); /* equality work vector */
294:     }
295:     if (tao->ineq_constrained) {
296:       VecDestroy(&auglag->Ps);     /* slack vars */
297:       PetscFree(auglag->Parr);     /* array of primal vectors */
298:       VecDestroy(&auglag->LgradS); /* slack grad */
299:       VecDestroy(&auglag->Cizero); /* zero vector for pointwise max */
300:       VecDestroy(&auglag->Yi);     /* inequality multipliers */
301:       VecDestroy(&auglag->Ciwork); /* inequality work vector */
302:       VecDestroy(&auglag->P);      /* combo primal solution */
303:       ISDestroy(&auglag->Pis[0]);  /* index set for X inside P */
304:       ISDestroy(&auglag->Pis[1]);  /* index set for S inside P */
305:       PetscFree(auglag->Pis);      /* array of P index sets */
306:       VecScatterDestroy(&auglag->Pscatter[0]);
307:       VecScatterDestroy(&auglag->Pscatter[1]);
308:       PetscFree(auglag->Pscatter);
309:       if (tao->eq_constrained) {
310:         VecDestroy(&auglag->Y);     /* combo multipliers */
311:         PetscFree(auglag->Yarr);    /* array of dual vectors */
312:         VecDestroy(&auglag->C);     /* combo constraints */
313:         ISDestroy(&auglag->Yis[0]); /* index set for Ye inside Y */
314:         ISDestroy(&auglag->Yis[1]); /* index set for Yi inside Y */
315:         PetscFree(auglag->Yis);
316:         VecScatterDestroy(&auglag->Yscatter[0]);
317:         VecScatterDestroy(&auglag->Yscatter[1]);
318:         PetscFree(auglag->Yscatter);
319:       }
320:     }
321:     if (tao->bounded) {
322:       VecDestroy(&auglag->PL); /* lower bounds for subsolver */
323:       VecDestroy(&auglag->PU); /* upper bounds for subsolver */
324:     }
325:   }
326:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", NULL);
327:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", NULL);
328:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", NULL);
329:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", NULL);
330:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", NULL);
331:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", NULL);
332:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", NULL);
333:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", NULL);
334:   PetscFree(tao->data);
335:   return 0;
336: }

338: static PetscErrorCode TaoSetFromOptions_ALMM(Tao tao, PetscOptionItems *PetscOptionsObject)
339: {
340:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
341:   PetscInt  i;

343:   PetscOptionsHeadBegin(PetscOptionsObject, "Augmented Lagrangian multiplier method solves problems with general constraints by converting them into a sequence of unconstrained problems.");
344:   PetscOptionsReal("-tao_almm_mu_init", "initial penalty parameter", "", auglag->mu0, &auglag->mu0, NULL);
345:   PetscOptionsReal("-tao_almm_mu_factor", "increase factor for the penalty parameter", "", auglag->mu_fac, &auglag->mu_fac, NULL);
346:   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:   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:   PetscOptionsReal("-tao_almm_mu_max", "maximum safeguard for penalty parameter updates", "", auglag->mu_max, &auglag->mu_max, NULL);
349:   PetscOptionsReal("-tao_almm_ye_min", "minimum safeguard for equality multiplier updates", "", auglag->ye_min, &auglag->ye_min, NULL);
350:   PetscOptionsReal("-tao_almm_ye_max", "maximum safeguard for equality multipliers updates", "", auglag->ye_max, &auglag->ye_max, NULL);
351:   PetscOptionsReal("-tao_almm_yi_min", "minimum safeguard for inequality multipliers updates", "", auglag->yi_min, &auglag->yi_min, NULL);
352:   PetscOptionsReal("-tao_almm_yi_max", "maximum safeguard for inequality multipliers updates", "", auglag->yi_max, &auglag->yi_max, NULL);
353:   PetscOptionsEnum("-tao_almm_type", "augmented Lagrangian formulation type for the subproblem", "TaoALMMType", TaoALMMTypes, (PetscEnum)auglag->type, (PetscEnum *)&auglag->type, NULL);
354:   PetscOptionsHeadEnd();
355:   TaoSetOptionsPrefix(auglag->subsolver, ((PetscObject)tao)->prefix);
356:   TaoAppendOptionsPrefix(auglag->subsolver, "tao_almm_subsolver_");
357:   TaoSetFromOptions(auglag->subsolver);
358:   for (i = 0; i < tao->numbermonitors; i++) {
359:     PetscObjectReference((PetscObject)tao->monitorcontext[i]);
360:     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:   return 0;
364: }

366: /* -------------------------------------------------------- */

368: /*MC
369:   TAOALMM - Augmented Lagrangian multiplier method for solving nonlinear optimization problems with general constraints.

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 <phr,classic>   - change formulation of the augmented Lagrangian merit function for the subproblem (default: phr)

383:   Level: beginner

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.

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 may converge
392:   faster for smaller problems but is highly susceptible to poor step lengths in the subproblem due to the positivity
393:   constraint on slack variables. PHR avoids this issue by eliminating the slack variables entirely, and is highly
394:   desirable for problems with a large number of inequality constraints.

396:   The subproblem is solved using a nested first-order TAO solver (default: `TAOBQNLS`). The user can retrieve a
397:   pointer to the subsolver via `TaoALMMGetSubsolver()` or pass command line arguments to it using the
398:   "-tao_almm_subsolver_" prefix. Currently, `TAOALMM` does not support second-order methods for the subproblem.

400: .vb
401:   while unconverged
402:     solve argmin_x L(x) s.t. l <= x <= u
403:     if ||c|| <= y_tol
404:       if ||c|| <= c_tol && ||Lgrad|| <= g_tol:
405:         problem converged, return solution
406:       else
407:         constraints sufficiently improved
408:         update multipliers and tighten tolerances
409:       endif
410:     else
411:       constraints did not improve
412:       update penalty and loosen tolerances
413:     endif
414:   endwhile
415: .ve

417: .seealso: `TAOALMM`, `Tao`, `TaoALMMGetType()`, `TaoALMMSetType()`, `TaoALMMSetSubsolver()`, `TaoALMMGetSubsolver()`,
418:           `TaoALMMGetMultipliers()`, `TaoALMMSetMultipliers()`, `TaoALMMGetPrimalIS()`, `TaoALMMGetDualIS()`
419: M*/
420: PETSC_EXTERN PetscErrorCode TaoCreate_ALMM(Tao tao)
421: {
422:   TAO_ALMM *auglag;

424:   PetscNew(&auglag);

426:   tao->ops->destroy        = TaoDestroy_ALMM;
427:   tao->ops->setup          = TaoSetUp_ALMM;
428:   tao->ops->setfromoptions = TaoSetFromOptions_ALMM;
429:   tao->ops->view           = TaoView_ALMM;
430:   tao->ops->solve          = TaoSolve_ALMM;

432:   tao->gatol = 1.e-5;
433:   tao->grtol = 0.0;
434:   tao->gttol = 0.0;
435:   tao->catol = 1.e-5;
436:   tao->crtol = 0.0;

438:   tao->data           = (void *)auglag;
439:   auglag->parent      = tao;
440:   auglag->mu0         = 10.0;
441:   auglag->mu          = auglag->mu0;
442:   auglag->mu_fac      = 10.0;
443:   auglag->mu_max      = PETSC_INFINITY;
444:   auglag->mu_pow_good = 0.9;
445:   auglag->mu_pow_bad  = 0.1;
446:   auglag->ye_min      = PETSC_NINFINITY;
447:   auglag->ye_max      = PETSC_INFINITY;
448:   auglag->yi_min      = PETSC_NINFINITY;
449:   auglag->yi_max      = PETSC_INFINITY;
450:   auglag->ytol0       = 0.1 / PetscPowReal(auglag->mu0, auglag->mu_pow_bad);
451:   auglag->ytol        = auglag->ytol0;
452:   auglag->gtol0       = 1.0 / auglag->mu0;
453:   auglag->gtol        = auglag->gtol0;

455:   auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
456:   auglag->type    = TAO_ALMM_PHR;
457:   auglag->info    = PETSC_FALSE;

459:   TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver);
460:   TaoSetType(auglag->subsolver, TAOBQNLS);
461:   TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0);
462:   TaoSetMaximumIterations(auglag->subsolver, 1000);
463:   TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000);
464:   TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY);
465:   PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1);

467:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private);
468:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private);
469:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private);
470:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private);
471:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private);
472:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private);
473:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private);
474:   PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private);
475:   return 0;
476: }

478: static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P)
479: {
480:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

482:   if (tao->ineq_constrained) {
483:     VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);
484:     VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);
485:     VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);
486:     VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);
487:   } else {
488:     VecCopy(X, P);
489:   }
490:   return 0;
491: }

493: static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y)
494: {
495:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

497:   if (tao->eq_constrained) {
498:     if (tao->ineq_constrained) {
499:       VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);
500:       VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);
501:       VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);
502:       VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);
503:     } else {
504:       VecCopy(EQ, Y);
505:     }
506:   } else {
507:     VecCopy(IN, Y);
508:   }
509:   return 0;
510: }

512: static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S)
513: {
514:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

516:   if (tao->ineq_constrained) {
517:     VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);
518:     VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);
519:     VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);
520:     VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);
521:   } else {
522:     VecCopy(P, X);
523:   }
524:   return 0;
525: }

527: /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */
528: static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao)
529: {
530:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

532:   /* if bounded, project the gradient */
533:   if (tao->bounded) VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX);
534:   if (auglag->type == TAO_ALMM_PHR) {
535:     VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm);
536:     auglag->cenorm = 0.0;
537:     if (tao->eq_constrained) VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm);
538:     auglag->cinorm = 0.0;
539:     if (tao->ineq_constrained) {
540:       VecCopy(auglag->Yi, auglag->Ciwork);
541:       VecScale(auglag->Ciwork, -1.0 / auglag->mu);
542:       VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork);
543:       VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm);
544:     }
545:     auglag->cnorm_old = auglag->cnorm;
546:     auglag->cnorm     = PetscMax(auglag->cenorm, auglag->cinorm);
547:     auglag->ytol      = auglag->ytol0 * auglag->cnorm_old;
548:   } else {
549:     VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm);
550:     VecNorm(auglag->C, NORM_2, &auglag->cnorm);
551:   }
552:   return 0;
553: }

555: static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao)
556: {
557:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

559:   /* split solution into primal and slack components */
560:   TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps);

562:   /* compute f, df/dx and the constraints */
563:   TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX);
564:   if (tao->eq_constrained) {
565:     TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce);
566:     TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae);
567:   }
568:   if (tao->ineq_constrained) {
569:     TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci);
570:     TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai);
571:     switch (auglag->type) {
572:     case TAO_ALMM_CLASSIC:
573:       /* classic formulation converts inequality to equality constraints via slack variables */
574:       VecAXPY(auglag->Ci, -1.0, auglag->Ps);
575:       break;
576:     case TAO_ALMM_PHR:
577:       /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */
578:       VecScale(auglag->Ci, -1.0);
579:       MatScale(auglag->Ai, -1.0);
580:       break;
581:     default:
582:       break;
583:     }
584:   }
585:   /* combine constraints into one vector */
586:   TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C);
587:   return 0;
588: }

590: /*
591: Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmin(0, Ci + Yi/mu)^T pmin(0, Ci + Yi/mu)]

593: dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmin(0, Ci + Yi/mu)^T Ai]

595: dLphr/dS = 0
596: */
597: static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao)
598: {
599:   TAO_ALMM *auglag  = (TAO_ALMM *)tao->data;
600:   PetscReal eq_norm = 0.0, ineq_norm = 0.0;

602:   TaoALMMEvaluateIterate_Private(tao);
603:   if (tao->eq_constrained) {
604:     /* Ce_work = mu*(Ce + Ye/mu) */
605:     VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce);
606:     VecDot(auglag->Cework, auglag->Cework, &eq_norm); /* contribution to scalar Lagrangian */
607:     VecScale(auglag->Cework, auglag->mu);
608:     /* dL/dX += mu*(Ce + Ye/mu)^T Ae */
609:     MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX);
610:   }
611:   if (tao->ineq_constrained) {
612:     /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */
613:     VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci);
614:     VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork);
615:     VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm); /* contribution to scalar Lagrangian */
616:     /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */
617:     VecScale(auglag->Ciwork, auglag->mu);
618:     MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX);
619:     /* dL/dS = 0 because there are no slacks in PHR */
620:     VecZeroEntries(auglag->LgradS);
621:   }
622:   /* combine gradient together */
623:   TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);
624:   /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */
625:   auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm);
626:   return 0;
627: }

629: /*
630: Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)]

632: dLc/dX = dF/dX + Ye^TAe + Yi^TAi + 0.5*mu*[Ce^TAe + (Ci - S)^TAi]

634: dLc/dS = -[Yi + mu*(Ci - S)]
635: */
636: static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao)
637: {
638:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
639:   PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0;

641:   TaoALMMEvaluateIterate_Private(tao);
642:   if (tao->eq_constrained) {
643:     /* compute scalar contributions */
644:     VecDot(auglag->Ye, auglag->Ce, &yeTce);
645:     VecDot(auglag->Ce, auglag->Ce, &ceTce);
646:     /* dL/dX += ye^T Ae */
647:     MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX);
648:     /* dL/dX += 0.5 * mu * ce^T Ae */
649:     MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork);
650:     VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork);
651:   }
652:   if (tao->ineq_constrained) {
653:     /* compute scalar contributions */
654:     VecDot(auglag->Yi, auglag->Ci, &yiTcims);
655:     VecDot(auglag->Ci, auglag->Ci, &cimsTcims);
656:     /* dL/dX += yi^T Ai */
657:     MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX);
658:     /* dL/dX += 0.5 * mu * (ci - s)^T Ai */
659:     MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork);
660:     VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork);
661:     /* dL/dS = -[yi + mu*(ci - s)] */
662:     VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi);
663:     VecScale(auglag->LgradS, -1.0);
664:   }
665:   /* combine gradient together */
666:   TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);
667:   /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */
668:   auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims);
669:   return 0;
670: }

672: PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, void *ctx)
673: {
674:   TAO_ALMM *auglag = (TAO_ALMM *)ctx;

676:   VecCopy(P, auglag->P);
677:   (*auglag->sub_obj)(auglag->parent);
678:   *Lval = auglag->Lval;
679:   return 0;
680: }

682: PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx)
683: {
684:   TAO_ALMM *auglag = (TAO_ALMM *)ctx;

686:   VecCopy(P, auglag->P);
687:   (*auglag->sub_obj)(auglag->parent);
688:   VecCopy(auglag->G, G);
689:   *Lval = auglag->Lval;
690:   return 0;
691: }