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, 0.0));
26: }
27: if (tao->eq_constrained) PetscCall(VecSet(auglag->Ye, 0.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: 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));
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: }
121: static PetscErrorCode TaoView_ALMM(Tao tao, PetscViewer viewer)
122: {
123: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
124: PetscBool isascii;
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: }
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;
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) PetscCall(VecDuplicate(auglag->Ci, &auglag->Cizero));
176: if (!auglag->Ps) { /* slack vars */
177: PetscCall(VecDuplicate(auglag->Ci, &auglag->Ps));
178: }
179: if (!auglag->LgradS) { /* slack component of Lagrangian gradient */
180: PetscCall(VecDuplicate(auglag->Ci, &auglag->LgradS));
181: }
182: /* create vector for combined primal space and the associated communication objects */
183: if (!auglag->P) {
184: PetscCall(PetscMalloc1(2, &auglag->Parr));
185: auglag->Parr[0] = auglag->Px;
186: auglag->Parr[1] = auglag->Ps;
187: PetscCall(VecConcatenate(2, auglag->Parr, &auglag->P, &auglag->Pis));
188: PetscCall(PetscMalloc1(2, &auglag->Pscatter));
189: PetscCall(VecScatterCreate(auglag->P, auglag->Pis[0], auglag->Px, NULL, &auglag->Pscatter[0]));
190: PetscCall(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: PetscCall(PetscMalloc1(2, &auglag->Yarr));
196: auglag->Yarr[0] = auglag->Ye;
197: auglag->Yarr[1] = auglag->Yi;
198: PetscCall(VecConcatenate(2, auglag->Yarr, &auglag->Y, &auglag->Yis));
199: PetscCall(PetscMalloc1(2, &auglag->Yscatter));
200: PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[0], auglag->Ye, NULL, &auglag->Yscatter[0]));
201: PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[1], auglag->Yi, NULL, &auglag->Yscatter[1]));
202: }
203: if (!auglag->C) PetscCall(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: PetscCall(VecDuplicate(auglag->P, &auglag->Psub));
216: PetscCall(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 != tao->default_gatol && tao->catol != tao->default_catol) {
224: PetscCall(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: PetscCall(TaoSetSolution(auglag->subsolver, auglag->Psub));
241: PetscCall(TaoSetObjective(auglag->subsolver, TaoALMMSubsolverObjective_Private, (void *)auglag));
242: PetscCall(TaoSetObjectiveAndGradient(auglag->subsolver, NULL, TaoALMMSubsolverObjectiveAndGradient_Private, (void *)auglag));
243: if (tao->bounded) {
244: /* make sure that the subsolver is a bound-constrained method */
245: PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg));
246: PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm));
247: if (is_cg) {
248: PetscCall(TaoSetType(auglag->subsolver, TAOBNCG));
249: PetscCall(PetscInfo(tao, "TAOCG detected for bound-constrained problem, switching to TAOBNCG.\n"));
250: }
251: if (is_lmvm) {
252: PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS));
253: PetscCall(PetscInfo(tao, "TAOLMVM detected for bound-constrained problem, switching to TAOBQNLS.\n"));
254: }
255: /* create lower and upper bound clone vectors for subsolver */
256: if (!auglag->PL) PetscCall(VecDuplicate(auglag->P, &auglag->PL));
257: if (!auglag->PU) PetscCall(VecDuplicate(auglag->P, &auglag->PU));
258: if (tao->ineq_constrained) {
259: /* create lower and upper bounds for slack, set lower to 0 */
260: PetscCall(VecDuplicate(auglag->Ci, &SL));
261: PetscCall(VecDuplicate(auglag->Ci, &SU));
262: PetscCall(VecSet(SU, PETSC_INFINITY));
263: /* combine opt var bounds with slack bounds */
264: PetscCall(TaoALMMCombinePrimal_Private(tao, tao->XL, SL, auglag->PL));
265: PetscCall(TaoALMMCombinePrimal_Private(tao, tao->XU, SU, auglag->PU));
266: /* destroy work vectors */
267: PetscCall(VecDestroy(&SL));
268: PetscCall(VecDestroy(&SU));
269: } else {
270: /* no inequality constraints, just copy bounds into the subsolver */
271: PetscCall(VecCopy(tao->XL, auglag->PL));
272: PetscCall(VecCopy(tao->XU, auglag->PU));
273: }
274: PetscCall(TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU));
275: } else {
276: /* CLASSIC's slack variable is bounded, so need to set bounds */
277: //TODO what happens for non-constrained ALMM CLASSIC?
278: if (auglag->type == TAO_ALMM_CLASSIC) {
279: if (tao->ineq_constrained) {
280: /* create lower and upper bound clone vectors for subsolver
281: * They should be NFINITY and INFINITY */
282: if (!auglag->PL) PetscCall(VecDuplicate(auglag->P, &auglag->PL));
283: if (!auglag->PU) PetscCall(VecDuplicate(auglag->P, &auglag->PU));
284: PetscCall(VecSet(auglag->PL, PETSC_NINFINITY));
285: PetscCall(VecSet(auglag->PU, PETSC_INFINITY));
286: /* create lower and upper bounds for slack, set lower to 0 */
287: PetscCall(VecDuplicate(auglag->Ci, &SL));
288: PetscCall(VecDuplicate(auglag->Ci, &SU));
289: PetscCall(VecSet(SU, PETSC_INFINITY));
290: /* PL, PU is already set. Only copy Slack variable parts */
291: PetscCall(VecScatterBegin(auglag->Pscatter[1], SL, auglag->PL, INSERT_VALUES, SCATTER_REVERSE));
292: PetscCall(VecScatterEnd(auglag->Pscatter[1], SL, auglag->PL, INSERT_VALUES, SCATTER_REVERSE));
293: PetscCall(VecScatterBegin(auglag->Pscatter[1], SU, auglag->PU, INSERT_VALUES, SCATTER_REVERSE));
294: PetscCall(VecScatterEnd(auglag->Pscatter[1], SU, auglag->PU, INSERT_VALUES, SCATTER_REVERSE));
295: /* destroy work vectors */
296: PetscCall(VecDestroy(&SL));
297: PetscCall(VecDestroy(&SU));
298: /* make sure that the subsolver is a bound-constrained method
299: * Unfortunately duplicate code */
300: PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg));
301: PetscCall(PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm));
302: if (is_cg) {
303: PetscCall(TaoSetType(auglag->subsolver, TAOBNCG));
304: PetscCall(PetscInfo(tao, "TAOCG detected for TAO_ALMM_CLASSIC, switching to TAOBNCG.\n"));
305: }
306: if (is_lmvm) {
307: PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS));
308: PetscCall(PetscInfo(tao, "TAOLMVM detected for TAO_ALMM_CLASSIC, switching to TAOBQNLS.\n"));
309: }
310: PetscCall(TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU));
311: }
312: }
313: }
314: PetscCall(TaoSetUp(auglag->subsolver));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode TaoDestroy_ALMM(Tao tao)
319: {
320: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
322: PetscFunctionBegin;
323: PetscCall(TaoDestroy(&auglag->subsolver));
324: PetscCall(VecDestroy(&auglag->Psub));
325: PetscCall(VecDestroy(&auglag->G));
326: if (tao->setupcalled) {
327: PetscCall(VecDestroy(&auglag->Xwork));
328: if (tao->eq_constrained) {
329: PetscCall(VecDestroy(&auglag->Ye)); /* equality multipliers */
330: PetscCall(VecDestroy(&auglag->Cework)); /* equality work vector */
331: }
332: if (tao->ineq_constrained) {
333: PetscCall(VecDestroy(&auglag->Ps)); /* slack vars */
334: PetscCall(PetscFree(auglag->Parr)); /* array of primal vectors */
335: PetscCall(VecDestroy(&auglag->LgradS)); /* slack grad */
336: PetscCall(VecDestroy(&auglag->Cizero)); /* zero vector for pointwise max */
337: PetscCall(VecDestroy(&auglag->Yi)); /* inequality multipliers */
338: PetscCall(VecDestroy(&auglag->Ciwork)); /* inequality work vector */
339: PetscCall(VecDestroy(&auglag->P)); /* combo primal solution */
340: PetscCall(ISDestroy(&auglag->Pis[0])); /* index set for X inside P */
341: PetscCall(ISDestroy(&auglag->Pis[1])); /* index set for S inside P */
342: PetscCall(PetscFree(auglag->Pis)); /* array of P index sets */
343: PetscCall(VecScatterDestroy(&auglag->Pscatter[0]));
344: PetscCall(VecScatterDestroy(&auglag->Pscatter[1]));
345: PetscCall(PetscFree(auglag->Pscatter));
346: if (tao->eq_constrained) {
347: PetscCall(VecDestroy(&auglag->Y)); /* combo multipliers */
348: PetscCall(PetscFree(auglag->Yarr)); /* array of dual vectors */
349: PetscCall(VecDestroy(&auglag->C)); /* combo constraints */
350: PetscCall(ISDestroy(&auglag->Yis[0])); /* index set for Ye inside Y */
351: PetscCall(ISDestroy(&auglag->Yis[1])); /* index set for Yi inside Y */
352: PetscCall(PetscFree(auglag->Yis));
353: PetscCall(VecScatterDestroy(&auglag->Yscatter[0]));
354: PetscCall(VecScatterDestroy(&auglag->Yscatter[1]));
355: PetscCall(PetscFree(auglag->Yscatter));
356: }
357: }
358: if (tao->bounded) {
359: PetscCall(VecDestroy(&auglag->PL)); /* lower bounds for subsolver */
360: PetscCall(VecDestroy(&auglag->PU)); /* upper bounds for subsolver */
361: } else {
362: if (auglag->type == TAO_ALMM_CLASSIC) {
363: PetscCall(VecDestroy(&auglag->PL)); /* lower bounds for subsolver */
364: PetscCall(VecDestroy(&auglag->PU)); /* upper bounds for subsolver */
365: }
366: }
367: }
368: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", NULL));
369: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", NULL));
370: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", NULL));
371: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", NULL));
372: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", NULL));
373: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", NULL));
374: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", NULL));
375: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", NULL));
376: PetscCall(PetscFree(tao->data));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode TaoSetFromOptions_ALMM(Tao tao, PetscOptionItems PetscOptionsObject)
381: {
382: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
384: PetscFunctionBegin;
385: PetscOptionsHeadBegin(PetscOptionsObject, "Augmented Lagrangian multiplier method solves problems with general constraints by converting them into a sequence of unconstrained problems.");
386: PetscCall(PetscOptionsReal("-tao_almm_mu_init", "initial penalty parameter", "", auglag->mu0, &auglag->mu0, NULL));
387: PetscCall(PetscOptionsReal("-tao_almm_mu_factor", "increase factor for the penalty parameter", "", auglag->mu_fac, &auglag->mu_fac, NULL));
388: 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));
389: 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));
390: PetscCall(PetscOptionsReal("-tao_almm_mu_max", "maximum safeguard for penalty parameter updates", "", auglag->mu_max, &auglag->mu_max, NULL));
391: PetscCall(PetscOptionsReal("-tao_almm_ye_min", "minimum safeguard for equality multiplier updates", "", auglag->ye_min, &auglag->ye_min, NULL));
392: PetscCall(PetscOptionsReal("-tao_almm_ye_max", "maximum safeguard for equality multipliers updates", "", auglag->ye_max, &auglag->ye_max, NULL));
393: PetscCall(PetscOptionsReal("-tao_almm_yi_min", "minimum safeguard for inequality multipliers updates", "", auglag->yi_min, &auglag->yi_min, NULL));
394: PetscCall(PetscOptionsReal("-tao_almm_yi_max", "maximum safeguard for inequality multipliers updates", "", auglag->yi_max, &auglag->yi_max, NULL));
395: PetscCall(PetscOptionsEnum("-tao_almm_type", "augmented Lagrangian formulation type for the subproblem", "TaoALMMType", TaoALMMTypes, (PetscEnum)auglag->type, (PetscEnum *)&auglag->type, NULL));
396: PetscOptionsHeadEnd();
397: PetscCall(TaoSetOptionsPrefix(auglag->subsolver, ((PetscObject)tao)->prefix));
398: PetscCall(TaoAppendOptionsPrefix(auglag->subsolver, "tao_almm_subsolver_"));
399: PetscCall(TaoSetFromOptions(auglag->subsolver));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*MC
404: TAOALMM - Augmented Lagrangian multiplier method for solving nonlinear optimization problems with general constraints.
406: Options Database Keys:
407: + -tao_almm_mu_init real - initial penalty parameter (default: 10.)
408: . -tao_almm_mu_factor real - increase factor for the penalty parameter (default: 100.)
409: . -tao_almm_mu_max real - maximum safeguard for penalty parameter updates (default: 1.e20)
410: . -tao_almm_mu_power_good real - exponential for penalty parameter when multiplier update is accepted (default: 0.9)
411: . -tao_almm_mu_power_bad real - exponential for penalty parameter when multiplier update is rejected (default: 0.1)
412: . -tao_almm_ye_min real - minimum safeguard for equality multiplier updates (default: -1.e20)
413: . -tao_almm_ye_max real - maximum safeguard for equality multiplier updates (default: 1.e20)
414: . -tao_almm_yi_min real - minimum safeguard for inequality multiplier updates (default: -1.e20)
415: . -tao_almm_yi_max real - maximum safeguard for inequality multiplier updates (default: 1.e20)
416: - -tao_almm_type (phr|classic) - change formulation of the augmented Lagrangian merit function for the subproblem (default: phr)
418: Level: beginner
420: Notes:
421: This method converts a constrained problem into a sequence of unconstrained problems via the augmented
422: Lagrangian merit function. Bound constraints are pushed down to the subproblem without any modifications.
424: Two formulations are offered for the subproblem: canonical Hestenes-Powell augmented Lagrangian with slack
425: variables for inequality constraints, and a slack-less Powell-Hestenes-Rockafellar (PHR) formulation utilizing a
426: pointwise max() penalty on inequality constraints. The canonical augmented Lagrangian formulation may converge
427: faster for smaller problems but is highly susceptible to poor step lengths in the subproblem due to the positivity
428: constraint on slack variables. PHR avoids this issue by eliminating the slack variables entirely, and is highly
429: desirable for problems with a large number of inequality constraints.
431: The subproblem is solved using a nested first-order TAO solver (default: `TAOBQNLS`). The user can retrieve a
432: pointer to the subsolver via `TaoALMMGetSubsolver()` or pass command line arguments to it using the
433: "-tao_almm_subsolver_" prefix. Currently, `TAOALMM` does not support second-order methods for the subproblem.
435: .vb
436: while unconverged
437: solve argmin_x L(x) s.t. l <= x <= u
438: if ||c|| <= y_tol
439: if ||c|| <= c_tol && ||Lgrad|| <= g_tol:
440: problem converged, return solution
441: else
442: constraints sufficiently improved
443: update multipliers and tighten tolerances
444: endif
445: else
446: constraints did not improve
447: update penalty and loosen tolerances
448: endif
449: endwhile
450: .ve
452: .seealso: `TAOALMM`, `Tao`, `TaoALMMGetType()`, `TaoALMMSetType()`, `TaoALMMSetSubsolver()`, `TaoALMMGetSubsolver()`,
453: `TaoALMMGetMultipliers()`, `TaoALMMSetMultipliers()`, `TaoALMMGetPrimalIS()`, `TaoALMMGetDualIS()`
454: M*/
455: PETSC_EXTERN PetscErrorCode TaoCreate_ALMM(Tao tao)
456: {
457: TAO_ALMM *auglag;
459: PetscFunctionBegin;
460: PetscCall(PetscNew(&auglag));
462: tao->ops->destroy = TaoDestroy_ALMM;
463: tao->ops->setup = TaoSetUp_ALMM;
464: tao->ops->setfromoptions = TaoSetFromOptions_ALMM;
465: tao->ops->view = TaoView_ALMM;
466: tao->ops->solve = TaoSolve_ALMM;
467: tao->uses_gradient = PETSC_TRUE;
469: PetscCall(TaoParametersInitialize(tao));
470: PetscObjectParameterSetDefault(tao, gatol, 1.e-5);
471: PetscObjectParameterSetDefault(tao, grtol, 0.0);
472: PetscObjectParameterSetDefault(tao, gttol, 0.0);
473: PetscObjectParameterSetDefault(tao, catol, 1.e-5);
474: PetscObjectParameterSetDefault(tao, crtol, 0.0);
476: tao->data = (void *)auglag;
477: auglag->parent = tao;
478: auglag->mu0 = 10.0;
479: auglag->mu = auglag->mu0;
480: auglag->mu_fac = 10.0;
481: auglag->mu_max = PETSC_INFINITY;
482: auglag->mu_pow_good = 0.9;
483: auglag->mu_pow_bad = 0.1;
484: auglag->ye_min = PETSC_NINFINITY;
485: auglag->ye_max = PETSC_INFINITY;
486: auglag->yi_min = PETSC_NINFINITY;
487: auglag->yi_max = PETSC_INFINITY;
488: auglag->ytol0 = 1.0 / PetscPowReal(auglag->mu0, auglag->mu_pow_bad);
489: auglag->ytol = auglag->ytol0;
490: auglag->gtol0 = 1.0 / auglag->mu0;
491: auglag->gtol = auglag->gtol0;
493: auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
494: auglag->type = TAO_ALMM_PHR;
496: PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver));
497: PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS));
498: PetscCall(TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0));
499: PetscCall(TaoSetMaximumIterations(auglag->subsolver, 1000));
500: PetscCall(TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000));
501: PetscCall(TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY));
502: PetscCall(PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1));
504: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private));
505: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private));
506: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private));
507: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private));
508: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private));
509: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private));
510: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private));
511: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P)
516: {
517: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
519: PetscFunctionBegin;
520: if (tao->ineq_constrained) {
521: PetscCall(VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE));
522: PetscCall(VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE));
523: PetscCall(VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE));
524: PetscCall(VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE));
525: } else {
526: PetscCall(VecCopy(X, P));
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y)
532: {
533: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
535: PetscFunctionBegin;
536: if (tao->eq_constrained) {
537: if (tao->ineq_constrained) {
538: PetscCall(VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE));
539: PetscCall(VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE));
540: PetscCall(VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE));
541: PetscCall(VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE));
542: } else {
543: PetscCall(VecCopy(EQ, Y));
544: }
545: } else {
546: PetscCall(VecCopy(IN, Y));
547: }
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S)
552: {
553: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
555: PetscFunctionBegin;
556: if (tao->ineq_constrained) {
557: PetscCall(VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD));
558: PetscCall(VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD));
559: PetscCall(VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD));
560: PetscCall(VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD));
561: } else {
562: PetscCall(VecCopy(P, X));
563: }
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */
568: static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao)
569: {
570: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
572: PetscFunctionBegin;
573: /* if bounded, project the gradient */
574: if (tao->bounded) PetscCall(VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX));
575: if (auglag->type == TAO_ALMM_PHR) {
576: PetscCall(VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm));
577: auglag->cenorm = 0.0;
578: if (tao->eq_constrained) PetscCall(VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm));
579: auglag->cinorm = 0.0;
580: if (tao->ineq_constrained) {
581: PetscCall(VecCopy(auglag->Yi, auglag->Ciwork));
582: PetscCall(VecScale(auglag->Ciwork, -1.0 / auglag->mu));
583: PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork));
584: PetscCall(VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm));
585: }
586: auglag->cnorm_old = auglag->cnorm;
587: auglag->cnorm = PetscMax(auglag->cenorm, auglag->cinorm);
588: auglag->ytol = auglag->ytol0 * auglag->cnorm_old;
589: } else {
590: PetscCall(VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm));
591: PetscCall(VecNorm(auglag->C, NORM_2, &auglag->cnorm));
592: }
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao)
597: {
598: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
600: PetscFunctionBegin;
601: /* split solution into primal and slack components */
602: PetscCall(TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps));
604: /* compute f, df/dx and the constraints */
605: PetscCall(TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX));
606: if (tao->eq_constrained) {
607: PetscCall(TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce));
608: PetscCall(TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae));
609: }
610: if (tao->ineq_constrained) {
611: PetscCall(TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci));
612: PetscCall(TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai));
613: switch (auglag->type) {
614: case TAO_ALMM_CLASSIC:
615: /* classic formulation converts inequality to equality constraints via slack variables */
616: PetscCall(VecAXPY(auglag->Ci, -1.0, auglag->Ps));
617: break;
618: case TAO_ALMM_PHR:
619: /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */
620: PetscCall(VecScale(auglag->Ci, -1.0));
621: PetscCall(MatScale(auglag->Ai, -1.0));
622: break;
623: default:
624: break;
625: }
626: }
627: /* combine constraints into one vector */
628: PetscCall(TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: /*
633: Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)]
635: dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmax(0, Ci + Yi/mu)^T Ai]
637: dLphr/dS = 0
638: */
639: static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao)
640: {
641: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
642: PetscReal eq_norm = 0.0, ineq_norm = 0.0;
644: PetscFunctionBegin;
645: PetscCall(TaoALMMEvaluateIterate_Private(tao));
646: if (tao->eq_constrained) {
647: /* Ce_work = mu*(Ce + Ye/mu) */
648: PetscCall(VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce));
649: PetscCall(VecDot(auglag->Cework, auglag->Cework, &eq_norm)); /* contribution to scalar Lagrangian */
650: PetscCall(VecScale(auglag->Cework, auglag->mu));
651: /* dL/dX += mu*(Ce + Ye/mu)^T Ae */
652: PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX));
653: }
654: if (tao->ineq_constrained) {
655: /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */
656: PetscCall(VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci));
657: PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork));
658: PetscCall(VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm)); /* contribution to scalar Lagrangian */
659: /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */
660: PetscCall(VecScale(auglag->Ciwork, auglag->mu));
661: PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX));
662: /* dL/dS = 0 because there are no slacks in PHR */
663: PetscCall(VecZeroEntries(auglag->LgradS));
664: }
665: /* combine gradient together */
666: PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G));
667: /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */
668: auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm);
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*
673: Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)]
675: dLc/dX = dF/dX + Ye^TAe + Yi^TAi + mu*[Ce^TAe + (Ci - S)^TAi]
677: dLc/dS = -[Yi + mu*(Ci - S)]
678: */
679: static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao)
680: {
681: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
682: PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0;
684: PetscFunctionBegin;
685: PetscCall(TaoALMMEvaluateIterate_Private(tao));
686: if (tao->eq_constrained) {
687: /* compute scalar contributions */
688: PetscCall(VecDot(auglag->Ye, auglag->Ce, &yeTce));
689: PetscCall(VecDot(auglag->Ce, auglag->Ce, &ceTce));
690: /* dL/dX += ye^T Ae */
691: PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX));
692: /* dL/dX += mu * ce^T Ae */
693: PetscCall(MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork));
694: PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork));
695: }
696: if (tao->ineq_constrained) {
697: /* compute scalar contributions */
698: PetscCall(VecDot(auglag->Yi, auglag->Ci, &yiTcims));
699: PetscCall(VecDot(auglag->Ci, auglag->Ci, &cimsTcims));
700: /* dL/dX += yi^T Ai */
701: PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX));
702: /* dL/dX += mu * (ci - s)^T Ai */
703: PetscCall(MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork));
704: PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork));
705: /* dL/dS = -[yi + mu*(ci - s)] */
706: PetscCall(VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi));
707: PetscCall(VecScale(auglag->LgradS, -1.0));
708: }
709: /* combine gradient together */
710: PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G));
711: /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */
712: auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims);
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, PetscCtx ctx)
717: {
718: TAO_ALMM *auglag = (TAO_ALMM *)ctx;
720: PetscFunctionBegin;
721: PetscCall(VecCopy(P, auglag->P));
722: PetscCall((*auglag->sub_obj)(auglag->parent));
723: *Lval = auglag->Lval;
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, PetscCtx ctx)
728: {
729: TAO_ALMM *auglag = (TAO_ALMM *)ctx;
731: PetscFunctionBegin;
732: PetscCall(VecCopy(P, auglag->P));
733: PetscCall((*auglag->sub_obj)(auglag->parent));
734: PetscCall(VecCopy(auglag->G, G));
735: *Lval = auglag->Lval;
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }