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