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 != tao->default_gatol && tao->catol != tao->default_catol) {
    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:   PetscCall(TaoParametersInitialize(tao));
    437:   PetscObjectParameterSetDefault(tao, gatol, 1.e-5);
    438:   PetscObjectParameterSetDefault(tao, grtol, 0.0);
    439:   PetscObjectParameterSetDefault(tao, gttol, 0.0);
    440:   PetscObjectParameterSetDefault(tao, catol, 1.e-5);
    441:   PetscObjectParameterSetDefault(tao, crtol, 0.0);

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

    460:   auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
    461:   auglag->type    = TAO_ALMM_PHR;
    462:   auglag->info    = PETSC_FALSE;

    464:   PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver));
    465:   PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS));
    466:   PetscCall(TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0));
    467:   PetscCall(TaoSetMaximumIterations(auglag->subsolver, 1000));
    468:   PetscCall(TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000));
    469:   PetscCall(TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY));
    470:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1));

    472:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private));
    473:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private));
    474:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private));
    475:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private));
    476:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private));
    477:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private));
    478:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private));
    479:   PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private));
    480:   PetscFunctionReturn(PETSC_SUCCESS);
    481: }

    483: static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P)
    484: {
    485:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

    487:   PetscFunctionBegin;
    488:   if (tao->ineq_constrained) {
    489:     PetscCall(VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE));
    490:     PetscCall(VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE));
    491:     PetscCall(VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE));
    492:     PetscCall(VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE));
    493:   } else {
    494:     PetscCall(VecCopy(X, P));
    495:   }
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }

    499: static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y)
    500: {
    501:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

    503:   PetscFunctionBegin;
    504:   if (tao->eq_constrained) {
    505:     if (tao->ineq_constrained) {
    506:       PetscCall(VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE));
    507:       PetscCall(VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE));
    508:       PetscCall(VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE));
    509:       PetscCall(VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE));
    510:     } else {
    511:       PetscCall(VecCopy(EQ, Y));
    512:     }
    513:   } else {
    514:     PetscCall(VecCopy(IN, Y));
    515:   }
    516:   PetscFunctionReturn(PETSC_SUCCESS);
    517: }

    519: static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S)
    520: {
    521:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

    523:   PetscFunctionBegin;
    524:   if (tao->ineq_constrained) {
    525:     PetscCall(VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD));
    526:     PetscCall(VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD));
    527:     PetscCall(VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD));
    528:     PetscCall(VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD));
    529:   } else {
    530:     PetscCall(VecCopy(P, X));
    531:   }
    532:   PetscFunctionReturn(PETSC_SUCCESS);
    533: }

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

    540:   PetscFunctionBegin;
    541:   /* if bounded, project the gradient */
    542:   if (tao->bounded) PetscCall(VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX));
    543:   if (auglag->type == TAO_ALMM_PHR) {
    544:     PetscCall(VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm));
    545:     auglag->cenorm = 0.0;
    546:     if (tao->eq_constrained) PetscCall(VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm));
    547:     auglag->cinorm = 0.0;
    548:     if (tao->ineq_constrained) {
    549:       PetscCall(VecCopy(auglag->Yi, auglag->Ciwork));
    550:       PetscCall(VecScale(auglag->Ciwork, -1.0 / auglag->mu));
    551:       PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork));
    552:       PetscCall(VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm));
    553:     }
    554:     auglag->cnorm_old = auglag->cnorm;
    555:     auglag->cnorm     = PetscMax(auglag->cenorm, auglag->cinorm);
    556:     auglag->ytol      = auglag->ytol0 * auglag->cnorm_old;
    557:   } else {
    558:     PetscCall(VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm));
    559:     PetscCall(VecNorm(auglag->C, NORM_2, &auglag->cnorm));
    560:   }
    561:   PetscFunctionReturn(PETSC_SUCCESS);
    562: }

    564: static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao)
    565: {
    566:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

    568:   PetscFunctionBegin;
    569:   /* split solution into primal and slack components */
    570:   PetscCall(TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps));

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

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

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

    605: dLphr/dS = 0
    606: */
    607: static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao)
    608: {
    609:   TAO_ALMM *auglag  = (TAO_ALMM *)tao->data;
    610:   PetscReal eq_norm = 0.0, ineq_norm = 0.0;

    612:   PetscFunctionBegin;
    613:   PetscCall(TaoALMMEvaluateIterate_Private(tao));
    614:   if (tao->eq_constrained) {
    615:     /* Ce_work = mu*(Ce + Ye/mu) */
    616:     PetscCall(VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce));
    617:     PetscCall(VecDot(auglag->Cework, auglag->Cework, &eq_norm)); /* contribution to scalar Lagrangian */
    618:     PetscCall(VecScale(auglag->Cework, auglag->mu));
    619:     /* dL/dX += mu*(Ce + Ye/mu)^T Ae */
    620:     PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX));
    621:   }
    622:   if (tao->ineq_constrained) {
    623:     /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */
    624:     PetscCall(VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci));
    625:     PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork));
    626:     PetscCall(VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm)); /* contribution to scalar Lagrangian */
    627:     /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */
    628:     PetscCall(VecScale(auglag->Ciwork, auglag->mu));
    629:     PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX));
    630:     /* dL/dS = 0 because there are no slacks in PHR */
    631:     PetscCall(VecZeroEntries(auglag->LgradS));
    632:   }
    633:   /* combine gradient together */
    634:   PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G));
    635:   /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */
    636:   auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm);
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

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

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

    645: dLc/dS = -[Yi + mu*(Ci - S)]
    646: */
    647: static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao)
    648: {
    649:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
    650:   PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0;

    652:   PetscFunctionBegin;
    653:   PetscCall(TaoALMMEvaluateIterate_Private(tao));
    654:   if (tao->eq_constrained) {
    655:     /* compute scalar contributions */
    656:     PetscCall(VecDot(auglag->Ye, auglag->Ce, &yeTce));
    657:     PetscCall(VecDot(auglag->Ce, auglag->Ce, &ceTce));
    658:     /* dL/dX += ye^T Ae */
    659:     PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX));
    660:     /* dL/dX += 0.5 * mu * ce^T Ae */
    661:     PetscCall(MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork));
    662:     PetscCall(VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork));
    663:   }
    664:   if (tao->ineq_constrained) {
    665:     /* compute scalar contributions */
    666:     PetscCall(VecDot(auglag->Yi, auglag->Ci, &yiTcims));
    667:     PetscCall(VecDot(auglag->Ci, auglag->Ci, &cimsTcims));
    668:     /* dL/dX += yi^T Ai */
    669:     PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX));
    670:     /* dL/dX += 0.5 * mu * (ci - s)^T Ai */
    671:     PetscCall(MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork));
    672:     PetscCall(VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork));
    673:     /* dL/dS = -[yi + mu*(ci - s)] */
    674:     PetscCall(VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi));
    675:     PetscCall(VecScale(auglag->LgradS, -1.0));
    676:   }
    677:   /* combine gradient together */
    678:   PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G));
    679:   /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */
    680:   auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims);
    681:   PetscFunctionReturn(PETSC_SUCCESS);
    682: }

    684: PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, void *ctx)
    685: {
    686:   TAO_ALMM *auglag = (TAO_ALMM *)ctx;

    688:   PetscFunctionBegin;
    689:   PetscCall(VecCopy(P, auglag->P));
    690:   PetscCall((*auglag->sub_obj)(auglag->parent));
    691:   *Lval = auglag->Lval;
    692:   PetscFunctionReturn(PETSC_SUCCESS);
    693: }

    695: PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx)
    696: {
    697:   TAO_ALMM *auglag = (TAO_ALMM *)ctx;

    699:   PetscFunctionBegin;
    700:   PetscCall(VecCopy(P, auglag->P));
    701:   PetscCall((*auglag->sub_obj)(auglag->parent));
    702:   PetscCall(VecCopy(auglag->G, G));
    703:   *Lval = auglag->Lval;
    704:   PetscFunctionReturn(PETSC_SUCCESS);
    705: }