MAST3
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
gcmma_interface.hpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #ifndef __mast_gcmma_optimization_interface_h__
21 #define __mast_gcmma_optimization_interface_h__
22 
23 // C++ includes
24 #include <iomanip>
25 
26 // MAST includes
27 #include <mast/base/mast_config.h>
29 #include <mast/base/exceptions.hpp>
31 
32 // libMesh includes
33 #include <libmesh/parallel.h>
34 
35 
36 extern "C" {
37 extern void raasta_(int *M, int *N,
38  double *RAA0, double *RAA,
39  double *XMIN, double *XMAX,
40  double *DF0DX, double *DFDX);
41 extern void asympg_(int *ITER, int *M, int *N,
42  double *ALBEFA, double *GHINIT,
43  double *GHDECR, double *GHINCR,
44  double *XVAL, double *XMIN, double *XMAX,
45  double *XOLD1, double *XOLD2,
46  double *XLOW, double *XUPP,
47  double *ALFA, double *BETA);
48 extern void mmasug_(int *ITER, int *M, int *N, double *GEPS, int *IYFREE,
49  double *XVAL, double *XMMA,
50  double *XMIN, double *XMAX,
51  double *XLOW, double *XUPP,
52  double *ALFA, double *BETA,
53  double *A, double *B, double *C, double *Y, double *Z,
54  double *RAA0, double *RAA, double *ULAM,
55  double *F0VAL, double *FVAL,
56  double *F0APP, double *FAPP,
57  double *FMAX, double *DF0DX, double *DFDX,
58  double *P, double *Q, double *P0, double *Q0,
59  double *UU, double *GRADF, double *DSRCH, double *HESSF);
60 extern void conser_(int *M, int *ICONSE,
61  double *GEPS, double *F0NEW, double *F0APP,
62  double *FNEW, double *FAPP);
63 extern void raaupd_(int *M, int *N, double *GEPS,
64  double *XMMA, double *XVAL,
65  double *XMIN, double *XMAX,
66  double *XLOW, double *XUPP,
67  double *F0NEW, double *FNEW,
68  double *F0APP, double *FAPP,
69  double *RAA0, double *RAA);
70 extern void xupdat_(int *N, int *ITER,
71  double *XMMA, double *XVAL,
72  double *XOLD1, double *XOLD2);
73 extern void fupdat_(int *M, double *F0NEW, double *FNEW,
74  double *F0VAL, double *FVAL);
75 }
76 
77 
78 namespace MAST {
79 namespace Optimization {
80 namespace Solvers {
81 
82 
83 
84 template <typename FunctionEvaluationType>
86 
87 public:
88 
89  GCMMAInterface(libMesh::Parallel::Communicator &comm):
90  _comm (comm),
91  constr_penalty (5.e1),
92  initial_rel_step (1.e-2),
93  asymptote_reduction (0.7),
94  asymptote_expansion (1.2),
95  rel_change_tol (1.e-8),
96  max_inner_iters (15),
97  max_iters (1000),
99  total_iter (0),
101  _feval (nullptr)
102  { }
103 
104 
105  virtual ~GCMMAInterface()
106  { }
107 
118 
119  inline void set_function_evaluation(FunctionEvaluationType& feval) {
120 
121  Assert0(!_feval, "Function evaluation object already set");
122 
123  _feval = &feval;
124  }
125 
126 
130  inline void init() {
131 
132  Assert0(_feval, "Function evaluation object not set");
133  int
134  N = _feval->n_vars(),
135  total_iter = 0;
136 
137  _XVAL.resize(N, 0.);
138  _XMIN.resize(N, 0.);
139  _XMAX.resize(N, 0.);
140 
141  _feval->init_dvar(_XVAL, _XMIN, _XMAX);
142  }
143 
144 
145  inline void optimize() {
146 
147 #if MAST_ENABLE_GCMMA == 1
148 
149  Assert0(_feval, "Function evaluation object not set");
150 
151  int
152  N = _feval->n_vars(),
153  M = _feval->n_eq() + _feval->n_ineq();
154 
155  Assert1(N > 0, N, "Design variables must be greater than 0");
156  Assert2(_XVAL.size() == N, _XVAL.size(), N, "Design variables must be initialized");
157  Assert2(_XMIN.size() == N, _XMIN.size(), N, "Design variables must be initialized");
158  Assert2(_XMAX.size() == N, _XMAX.size(), N, "Design variables must be initialized");
159 
160  std::vector<real_t> XOLD1(N, 0.), XOLD2(N, 0.),
161  XMMA(N, 0.), XLOW(N, 0.), XUPP(N, 0.),
162  ALFA(N, 0.), BETA(N, 0.), DF0DX(N, 0.),
163  A(M, 0.), B(M, 0.), C(M, 0.), Y(M, 0.), RAA(M, 0.), ULAM(M, 0.),
164  FVAL(M, 0.), FAPP(M, 0.), FNEW(M, 0.), FMAX(M, 0.),
165  DFDX(M*N, 0.), P(M*N, 0.), Q(M*N, 0.), P0(N, 0.), Q0(N, 0.),
166  UU(M, 0.), GRADF(M, 0.), DSRCH(M, 0.), HESSF(M*(M+1)/2, 0.),
167  f0_iters(n_rel_change_iters);
168 
169  std::vector<int> IYFREE(M, 0);
170  std::vector<bool> eval_grads(M, false);
171 
172  real_t
173  ALBEFA = 0.1,
174  GHINIT = initial_rel_step,
175  GHDECR = asymptote_reduction,
176  GHINCR = asymptote_expansion,
177  F0VAL = 0.,
178  F0NEW = 0.,
179  F0APP = 0.,
180  RAA0 = 0.,
181  Z = 0.,
182  GEPS = rel_change_tol;
183 
184 
185  /*C********+*********+*********+*********+*********+*********+*********+
186  C
187  C The meaning of some of the scalars and vectors in the program:
188  C
189  C N = Complex of variables x_j in the problem.
190  C M = Complex of constraints in the problem (not including
191  C the simple upper and lower bounds on the variables).
192  C ALBEFA = Relative spacing between asymptote and mode limit. Lower value
193  C will cause the move limit (alpha,beta) to move closer to asymptote
194  C values (l, u).
195  C GHINIT = Initial asymptote setting. For the first two iterations the
196  C asymptotes (l, u) are defined based on offsets from the design
197  C point as this fraction of the design variable bounds, ie.
198  C l_j = x_j^k - GHINIT * (x_j^max - x_j^min)
199  C u_j = x_j^k + GHINIT * (x_j^max - x_j^min)
200  C GHDECR = Fraction by which the asymptote is reduced for oscillating
201  C changes in design variables based on three consecutive iterations
202  C GHINCR = Fraction by which the asymptote is increased for non-oscillating
203  C changes in design variables based on three consecutive iterations
204  C INNMAX = Maximal number of inner iterations within each outer iter.
205  C A reasonable choice is INNMAX=10.
206  C ITER = Current outer iteration number ( =1 the first iteration).
207  C GEPS = Tolerance parameter for the constraints.
208  C (Used in the termination criteria for the subproblem.)
209  C
210  C XVAL(j) = Current value of the variable x_j.
211  C XOLD1(j) = Value of the variable x_j one iteration ago.
212  C XOLD2(j) = Value of the variable x_j two iterations ago.
213  C XMMA(j) = Optimal value of x_j in the MMA subproblem.
214  C XMIN(j) = Original lower bound for the variable x_j.
215  C XMAX(j) = Original upper bound for the variable x_j.
216  C XLOW(j) = Value of the lower asymptot l_j.
217  C XUPP(j) = Value of the upper asymptot u_j.
218  C ALFA(j) = Lower bound for x_j in the MMA subproblem.
219  C BETA(j) = Upper bound for x_j in the MMA subproblem.
220  C F0VAL = Value of the objective function f_0(x)
221  C FVAL(i) = Value of the i:th constraint function f_i(x).
222  C DF0DX(j) = Derivative of f_0(x) with respect to x_j.
223  C FMAX(i) = Right hand side of the i:th constraint.
224  C DFDX(k) = Derivative of f_i(x) with respect to x_j,
225  C where k = (j-1)*M + i.
226  C P(k) = Coefficient p_ij in the MMA subproblem, where
227  C k = (j-1)*M + i.
228  C Q(k) = Coefficient q_ij in the MMA subproblem, where
229  C k = (j-1)*M + i.
230  C P0(j) = Coefficient p_0j in the MMA subproblem.
231  C Q0(j) = Coefficient q_0j in the MMA subproblem.
232  C B(i) = Right hand side b_i in the MMA subproblem.
233  C F0APP = Value of the approximating objective function
234  C at the optimal soultion of the MMA subproblem.
235  C FAPP(i) = Value of the approximating i:th constraint function
236  C at the optimal soultion of the MMA subproblem.
237  C RAA0 = Parameter raa_0 in the MMA subproblem.
238  C RAA(i) = Parameter raa_i in the MMA subproblem.
239  C Y(i) = Value of the "artificial" variable y_i.
240  C Z = Value of the "minimax" variable z.
241  C A(i) = Coefficient a_i for the variable z.
242  C C(i) = Coefficient c_i for the variable y_i.
243  C ULAM(i) = Value of the dual variable lambda_i.
244  C GRADF(i) = Gradient component of the dual objective function.
245  C DSRCH(i) = Search direction component in the dual subproblem.
246  C HESSF(k) = Hessian matrix component of the dual function.
247  C IYFREE(i) = 0 for dual variables which are fixed to zero in
248  C the current subspace of the dual subproblem,
249  C = 1 for dual variables which are "free" in
250  C the current subspace of the dual subproblem.
251  C
252  C********+*********+*********+*********+*********+*********+*********+*/
253 
254 
255  /*
256  * The USER should now give values to the parameters
257  * M, N, GEPS, XVAL (starting point),
258  * XMIN, XMAX, FMAX, A and C.
259  */
260  // _initi(M,N,GEPS,XVAL,XMIN,XMAX,FMAX,A,C);
261  // Assumed: FMAX == A
262  //_feval->init_dvar(XVAL, XMIN, XMAX);
263  // set the value of C[i] to be very large numbers
264  real_t max_x = 0.;
265  for (uint_t i=0; i<N; i++)
266  if (max_x < fabs(_XVAL[i]))
267  max_x = fabs(_XVAL[i]);
268 
269  int INNMAX=max_inner_iters, ITER=0, ITE=0, INNER=0, ICONSE=0;
270  /*
271  * The outer iterative process starts.
272  */
273  bool terminate = false, inner_terminate=false;
274  while (!terminate) {
275 
276  std::fill(C.begin(), C.end(), std::max(1.e0*max_x, constr_penalty));
277  GHINIT = initial_rel_step,
278  GHDECR = asymptote_reduction,
279  GHINCR = asymptote_expansion,
280 
281  total_iter++;
282  ITER=ITER+1;
283  ITE=ITE+1;
284  /*
285  * The USER should now calculate function values and gradients
286  * at XVAL. The result should be put in F0VAL,DF0DX,FVAL,DFDX.
287  */
288  std::fill(eval_grads.begin(), eval_grads.end(), true);
290  F0VAL, true, DF0DX,
291  FVAL, eval_grads, DFDX);
292  if (ITER == 1)
293  // output the very first iteration
294  _feval->output(total_iter, _XVAL, F0VAL, FVAL);
295 
296  /*
297  * RAA0,RAA,XLOW,XUPP,ALFA and BETA are calculated.
298  */
299  raasta_(&M, &N, &RAA0, &RAA[0], &_XMIN[0], &_XMAX[0], &DF0DX[0], &DFDX[0]);
300  asympg_(&ITER, &M, &N, &ALBEFA, &GHINIT, &GHDECR, &GHINCR,
301  &_XVAL[0], &_XMIN[0], &_XMAX[0], &XOLD1[0], &XOLD2[0],
302  &XLOW[0], &XUPP[0], &ALFA[0], &BETA[0]);
303  /*
304  * The inner iterative process starts.
305  */
306 
307  // write the asymptote data for the inneriterations
308  if (write_internal_iteration_data)
309  _output_iteration_data(ITER, _XVAL, _XMIN, _XMAX, XLOW, XUPP, ALFA, BETA);
310 
311  INNER=0;
312  inner_terminate = false;
313  while (!inner_terminate) {
314 
315  /*
316  * The subproblem is generated and solved.
317  */
318  mmasug_(&ITER, &M, &N, &GEPS, &IYFREE[0], &_XVAL[0], &XMMA[0],
319  &_XMIN[0], &_XMAX[0], &XLOW[0], &XUPP[0], &ALFA[0], &BETA[0],
320  &A[0], &B[0], &C[0], &Y[0], &Z, &RAA0, &RAA[0], &ULAM[0],
321  &F0VAL, &FVAL[0], &F0APP, &FAPP[0], &FMAX[0], &DF0DX[0], &DFDX[0],
322  &P[0], &Q[0], &P0[0], &Q0[0], &UU[0], &GRADF[0], &DSRCH[0], &HESSF[0]);
323  /*
324  * The USER should now calculate function values at XMMA.
325  * The result should be put in F0NEW and FNEW.
326  */
327  std::fill(eval_grads.begin(), eval_grads.end(), false);
328  _evaluate_wrapper(XMMA,
329  F0NEW, false, DF0DX,
330  FNEW, eval_grads, DFDX);
331 
333  // if the solution is poor, backtrack
334  std::vector<real_t> XMMA_new(XMMA);
335  real_t frac = 0.5;
336  while (M && FNEW[0] > 1.e2) {
337  std::cout << "*** Backtracking: frac = "
338  << frac
339  << " constr: " << FNEW[0]
340  << std::endl;
341  for (uint_t i=0; i<XMMA.size(); i++)
342  XMMA_new[i] = XOLD1[i] + frac*(XMMA[i]-XOLD1[i]);
343 
344  _evaluate_wrapper(XMMA_new,
345  F0NEW, false, DF0DX,
346  FNEW, eval_grads, DFDX);
347  frac *= frac;
348  }
349  for (uint_t i=0; i<XMMA.size(); i++)
350  XMMA[i] = XMMA_new[i];
352 
353  if (INNER >= INNMAX) {
354  std::cout
355  << "** Max Inner Iter Reached: Terminating! Inner Iter = "
356  << INNER << std::endl;
357  inner_terminate = true;
358  }
359  else {
360  /*
361  * It is checked if the approximations were conservative.
362  */
363  conser_( &M, &ICONSE, &GEPS, &F0NEW, &F0APP, &FNEW[0], &FAPP[0]);
364  if (ICONSE == 1) {
365  std::cout
366  << "** Conservative Solution: Terminating! Inner Iter = "
367  << INNER << std::endl;
368  inner_terminate = true;
369  }
370  else {
371  /*
372  * The approximations were not conservative, so RAA0 and RAA
373  * are updated and one more inner iteration is started.
374  */
375  INNER=INNER+1;
376  raaupd_( &M, &N, &GEPS, &XMMA[0], &_XVAL[0],
377  &_XMIN[0], &_XMAX[0], &XLOW[0], &XUPP[0],
378  &F0NEW, &FNEW[0], &F0APP, &FAPP[0], &RAA0, &RAA[0]);
379  }
380  }
381  }
382 
383  /*
384  * The inner iterative process has terminated, which means
385  * that an outer iteration has been completed.
386  * The variables are updated so that XVAL stands for the new
387  * outer iteration point. The fuction values are also updated.
388  */
389  xupdat_( &N, &ITER, &XMMA[0], &_XVAL[0], &XOLD1[0], &XOLD2[0]);
390  fupdat_( &M, &F0NEW, &FNEW[0], &F0VAL, &FVAL[0]);
391  /*
392  * The USER may now write the current solution.
393  */
394  _feval->output(total_iter, _XVAL, F0VAL, FVAL);
395  f0_iters[(ITE-1)%n_rel_change_iters] = F0VAL;
396 
397  /*
398  * One more outer iteration is started as long as
399  * ITE is less than MAXITE:
400  */
401  if (ITE == max_iters) {
402  std::cout
403  << "GCMMA: Reached maximum iterations, terminating! "
404  << std::endl;
405  terminate = true;
406  }
407 
408  // relative change in objective
409  bool rel_change_conv = true;
410  real_t f0_curr = f0_iters[n_rel_change_iters-1];
411 
412  for (uint_t i=0; i<n_rel_change_iters-1; i++) {
413  if (f0_curr > sqrt(GEPS))
414  rel_change_conv = (rel_change_conv &&
415  fabs(f0_iters[i]-f0_curr)/fabs(f0_curr) < GEPS);
416  else
417  rel_change_conv = (rel_change_conv &&
418  fabs(f0_iters[i]-f0_curr) < GEPS);
419  }
420  if (rel_change_conv) {
421  std::cout
422  << "GCMMA: Converged relative change tolerance, terminating! "
423  << std::endl;
424  terminate = true;
425  }
426  // tell all processors about the decision here
427  _comm.broadcast(terminate, 0);
428  }
429 
430 #endif //MAST_ENABLE_GCMMA == 1
431  }
432 
433 
434 private:
435 
436  libMesh::Parallel::Communicator &_comm;
437 
439  const std::vector<real_t>& XVAL,
440  const std::vector<real_t>& XMIN,
441  const std::vector<real_t>& XMAX,
442  const std::vector<real_t>& XLOW,
443  const std::vector<real_t>& XUPP,
444  const std::vector<real_t>& ALFA,
445  const std::vector<real_t>& BETA) {
446 
447  Assert0(_feval, "Evaluation function not initialized");
448  Assert2(XVAL.size() == _feval->n_vars(),
449  XVAL.size(), _feval->n_vars(),
450  "Incorrect vector size");
451  Assert2(XMIN.size() == _feval->n_vars(),
452  XMIN.size(), _feval->n_vars(),
453  "Incorrect vector size");
454  Assert2(XMAX.size() == _feval->n_vars(),
455  XMAX.size(), _feval->n_vars(),
456  "Incorrect vector size");
457  Assert2(XLOW.size() == _feval->n_vars(),
458  XLOW.size(), _feval->n_vars(),
459  "Incorrect vector size");
460  Assert2(XUPP.size() == _feval->n_vars(),
461  XUPP.size(), _feval->n_vars(),
462  "Incorrect vector size");
463  Assert2(ALFA.size() == _feval->n_vars(),
464  ALFA.size(), _feval->n_vars(),
465  "Incorrect vector size");
466  Assert2(BETA.size() == _feval->n_vars(),
467  BETA.size(), _feval->n_vars(),
468  "Incorrect vector size");
469 
470  std::cout
471  << "****************************************************\n"
472  << " GCMMA: ASYMPTOTE DATA \n"
473  << "****************************************************\n"
474  << std::setw(5) << "Iter: " << i << std::endl
475  << std::setw(5) << "DV"
476  << std::setw(20) << "XMIN"
477  << std::setw(20) << "XLOW"
478  << std::setw(20) << "ALFA"
479  << std::setw(20) << "X"
480  << std::setw(20) << "BETA"
481  << std::setw(20) << "XUP"
482  << std::setw(20) << "XMAX" << std::endl;
483 
484  for (uint_t j=0; j<_feval->n_vars(); j++)
485  std::cout
486  << std::setw(5) << j
487  << std::setw(20) << XMIN[j]
488  << std::setw(20) << XLOW[j]
489  << std::setw(20) << ALFA[j]
490  << std::setw(20) << XVAL[j]
491  << std::setw(20) << BETA[j]
492  << std::setw(20) << XUPP[j]
493  << std::setw(20) << XMAX[j] << std::endl;
494  }
495 
496  inline void
497  _evaluate_wrapper(std::vector<real_t> &x,
498  real_t &obj,
499  bool eval_obj_grad,
500  std::vector<real_t> &obj_grad,
501  std::vector<real_t> &fvals,
502  std::vector<bool> &eval_grads,
503  std::vector<real_t> &grads) {
504 
505  // rank 0 will broadcase the DV values to all ranks
506  _comm.broadcast(x, 0/*, true*/);
507  _comm.broadcast(eval_obj_grad, 0);
508  //_comm.broadcast(eval_grads, 0/*, true*/);
509 
510  _feval->evaluate(x,
511  obj,
512  eval_obj_grad,
513  obj_grad,
514  fvals,
515  eval_grads,
516  grads);
517 
518  // make sure that all data returned is consistent across
519  // processors
520  Assert0(_comm.verify(obj),
521  "Objective function has different values on ranks");
522  Assert0(_comm.verify(obj_grad),
523  "Objective function gradient has different values on ranks");
524  Assert0(_comm.verify(fvals),
525  "Constraint functions has different values on ranks");
526  Assert0(_comm.verify(grads),
527  "Constraint function gradient has different values on ranks");
528  }
529 
530  FunctionEvaluationType *_feval;
531  std::vector<real_t> _XVAL, _XMIN, _XMAX;
532 };
533 
534 } // namespace Solvers
535 } // namespace Optimization
536 } // namespace MAST
537 
538 #endif // __mast_gcmma_optimization_interface_h__
void _evaluate_wrapper(std::vector< real_t > &x, real_t &obj, bool eval_obj_grad, std::vector< real_t > &obj_grad, std::vector< real_t > &fvals, std::vector< bool > &eval_grads, std::vector< real_t > &grads)
void raasta_(int *M, int *N, double *RAA0, double *RAA, double *XMIN, double *XMAX, double *DF0DX, double *DFDX)
void fupdat_(int *M, double *F0NEW, double *FNEW, double *F0VAL, double *FVAL)
void mmasug_(int *ITER, int *M, int *N, double *GEPS, int *IYFREE, double *XVAL, double *XMMA, double *XMIN, double *XMAX, double *XLOW, double *XUPP, double *ALFA, double *BETA, double *A, double *B, double *C, double *Y, double *Z, double *RAA0, double *RAA, double *ULAM, double *F0VAL, double *FVAL, double *F0APP, double *FAPP, double *FMAX, double *DF0DX, double *DFDX, double *P, double *Q, double *P0, double *Q0, double *UU, double *GRADF, double *DSRCH, double *HESSF)
libMesh::Parallel::Communicator & _comm
#define Assert1(cond, v1, msg)
Definition: exceptions.hpp:143
void conser_(int *M, int *ICONSE, double *GEPS, double *F0NEW, double *F0APP, double *FNEW, double *FAPP)
void xupdat_(int *N, int *ITER, double *XMMA, double *XVAL, double *XOLD1, double *XOLD2)
void init()
initializes the design variable vector from the function object.
void raaupd_(int *M, int *N, double *GEPS, double *XMMA, double *XVAL, double *XMIN, double *XMAX, double *XLOW, double *XUPP, double *F0NEW, double *FNEW, double *F0APP, double *FAPP, double *RAA0, double *RAA)
#define Assert0(cond, msg)
Definition: exceptions.hpp:134
unsigned int uint_t
#define Assert2(cond, v1, v2, msg)
Definition: exceptions.hpp:152
double real_t
void asympg_(int *ITER, int *M, int *N, double *ALBEFA, double *GHINIT, double *GHDECR, double *GHINCR, double *XVAL, double *XMIN, double *XMAX, double *XOLD1, double *XOLD2, double *XLOW, double *XUPP, double *ALFA, double *BETA)
void set_function_evaluation(FunctionEvaluationType &feval)
void _output_iteration_data(uint_t i, const std::vector< real_t > &XVAL, const std::vector< real_t > &XMIN, const std::vector< real_t > &XMAX, const std::vector< real_t > &XLOW, const std::vector< real_t > &XUPP, const std::vector< real_t > &ALFA, const std::vector< real_t > &BETA)
GCMMAInterface(libMesh::Parallel::Communicator &comm)