/* nag_glopt_nlp_multistart_sqp (e05ucc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void schwefel_obj(Integer *mode, Integer n, const double *x, double *objf, double *objgrd, Integer nstate, Nag_Comm *comm); static void schwefel_confun(Integer *mode, Integer ncnln, Integer n, Integer tdcjsl, const Integer *needc, const double *x, double *c, double *cjsl, Integer nstate, Nag_Comm *comm); static void mystart(Integer npts, double quas[], Integer n, Nag_Boolean repeat, const double bl[], const double bu[], Nag_Comm *comm, Integer *mode); #ifdef __cplusplus } #endif int main(void) { Integer exit_status = 0; Integer print_all_solutions = 0; Integer liopts = 740, lopts = 485, n = 2, nclin = 1, ncnln = 2; /* Scalars */ Integer i, ic, j, l, nb, npts, tda, ldcjac, sdcjac, ldr, sdr, ldx, ldobjgrd, ldclamda, ldistate, ldc; /* Arrays */ double *a=0, *bl=0, *bu=0, *c=0, *cjac=0, *clamda=0, *objf=0, *objgrd=0, *r=0, *opts=0, *work=0, *x=0; Integer *info=0, *istate=0, *iter=0, *iopts=0; char nag_enum_arg[40]; /* NAG Types */ NagError fail; Nag_Comm comm; Nag_Boolean repeat; INIT_FAIL(fail); printf("nag_glopt_nlp_multistart_sqp (e05ucc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%ld%ld%*[^\n]", &nb, &npts); scanf("%s%*[^\n]", nag_enum_arg); /* nag_enum_name_to_value (x04nac). * Converts NAG enum member name to value */ repeat = (Nag_Boolean) nag_enum_name_to_value(nag_enum_arg); /* The minimum trailing dimension for a is tda = n (or 1). */ if (nclin>0) { tda = n; if (!(a = NAG_ALLOC(nclin*tda, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else tda = 1; #define A(I,J) a[(I-1)*tda + (J-1)] #define X(I,J) x[(J-1)*ldx + (I-1)] #define ISTATE(I,J) istate[(J-1)*ldistate + (I-1)] #define CLAMDA(I,J) clamda[(J-1)*ldclamda + (I-1)] #define C(I,J) c[(J-1)*ldc + (I-1)] ldx = n; ldobjgrd = n; ldc = ncnln; ldcjac = ncnln; if (ncnln>0) { sdcjac = n; if ( !(c = NAG_ALLOC(ldc*nb, double)) || !(cjac = NAG_ALLOC(ldcjac*sdcjac*nb, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else sdcjac = 0; ldr = n; sdr = n; ldclamda = n + nclin + ncnln; ldistate = n + nclin + ncnln; if ( !(bl = NAG_ALLOC(n + nclin + ncnln, double)) || !(bu = NAG_ALLOC(n + nclin + ncnln, double)) || !(clamda = NAG_ALLOC(ldclamda*nb, double)) || !(objf = NAG_ALLOC(nb, double)) || !(objgrd = NAG_ALLOC(ldobjgrd*nb, double)) || !(r = NAG_ALLOC(ldr*sdr*nb, double)) || !(opts = NAG_ALLOC(lopts, double)) || !(work = NAG_ALLOC(nclin, double)) || !(x = NAG_ALLOC(ldx*nb, double)) || !(info = NAG_ALLOC(nb, Integer)) || !(istate = NAG_ALLOC(ldistate*nb, Integer)) || !(iter = NAG_ALLOC(nb, Integer)) || !(iopts = NAG_ALLOC(liopts, Integer)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } bl[0] = -500.0; bl[1] = -500.0; bl[2] = -10000.0; bl[3] = -1.0; bl[4] = -0.9; bu[0] = 500.0; bu[1] = 500.0; bu[2] = 10.0; bu[3] = 500000.0; bu[4] = 0.9; A(1, 1) = 3.0; A(1, 2) = -2.0; /* Initialize nag_glopt_nlp_multistart_sqp (e05ucc). * nag_glopt_opt_set (e05zkc). * Option setting routine for global optimization. */ nag_glopt_opt_set("Initialize = e05ucc", iopts, liopts, opts, lopts, NAGERR_DEFAULT); /* Solve the problem with repeatable random starting points using * nag_glopt_nlp_multistart_sqp (e05ucc). * Global optimization using multi-start, nonlinear constraints. * Pass 'mystart' instead of NULLFN as the start procedure to set * your own start points. */ nag_glopt_nlp_multistart_sqp(n, nclin, ncnln, a, tda, bl, bu, schwefel_confun, schwefel_obj, npts, x, ldx, NULLFN, repeat, nb, objf, objgrd, ldobjgrd, iter, c, ldc, cjac, ldcjac, sdcjac, r, ldr, sdr, clamda, ldclamda, istate, ldistate, iopts, opts, &comm, info, &fail); /* Check for error exits. */ switch (fail.code) { case NE_NOERROR: l = nb; break; case NW_SOME_SOLUTIONS: l = info[nb-1]; printf("Only %ld solutions found\n", l); break; default: exit_status = 2; printf("Error from nag_glopt_nlp_multistart_sqp (e05ucc)\n%s\n", fail.message); goto END; } for (i=1; i<=l; i++) { printf("Solution number %ld\n\n", i); printf("Local minimization exited with code %ld\n", info[i-1]); printf("\nVarbl Istate Value Lagr Mult\n\n\n"); for (j=1; j<=n; j++) printf("V %3ld %3ld %14.6g %12.4g\n", j, ISTATE(j,i), X(j,i), CLAMDA(j,i)); if (nclin>0) { printf("\nL Con Istate Value Lagr Mult\n\n"); /* nag_dgemv (f16pac) performs the matrix vector multiplication A*x * (linear constraint values) and puts the result in * the first nclin locations of work. */ nag_dgemv(Nag_RowMajor, Nag_NoTrans, nclin, n, 1.0, a, tda, &X(1,i), 1, 0.0, work, 1, NAGERR_DEFAULT); for (j = n+1; j <= n+nclin; j++) printf("L %3ld %3ld %14.6g %12.4g\n", j-n, ISTATE(j,i), work[j-n-1], CLAMDA(j,i)); } if (ncnln>0) { printf("\n\nN Con Istate Value Lagr Mult\n\n"); for (j = n+nclin+1; j <= n+nclin+ncnln; j++) { ic = j - n - nclin; printf("N %3ld %3ld %14.6g %12.4g\n", ic, ISTATE(j,i), C(ic,i), CLAMDA(j,i)); } } printf("\n\nFinal objective value = %15.7g\n", objf[i-1]); printf("\nQP multipliers\n"); for (j = 1; j <= n+nclin+ncnln; j++) printf("%12.4e\n", CLAMDA(j,i)); if (l==1) goto END; if (print_all_solutions==0) { printf("\n(Printing of further solutions suppressed)\n"); goto END; } printf("\n"); for (j = 0; j < 61; j++) printf("-"); printf("\n"); } END: if (a) NAG_FREE(a); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); if (c) NAG_FREE(c); if (cjac) NAG_FREE(cjac); if (clamda) NAG_FREE(clamda); if (objf) NAG_FREE(objf); if (objgrd) NAG_FREE(objgrd); if (r) NAG_FREE(r); if (opts) NAG_FREE(opts); if (work) NAG_FREE(work); if (x) NAG_FREE(x); if (info) NAG_FREE(info); if (istate) NAG_FREE(istate); if (iter) NAG_FREE(iter); if (iopts) NAG_FREE(iopts); return exit_status; } static void schwefel_obj(Integer *mode, Integer n, const double *x, double *objf, double *objgrd, Integer nstate, Nag_Comm *comm) { /* Scalars */ double t; if (*mode==0 || *mode==2) { /* Evaluate the objective function. */ *objf = x[0]*sin(sqrt(fabs(x[0]))) + x[1]*sin(sqrt(fabs(x[1]))); } if (*mode==1 || *mode==2) { /* Calculate the gradient of the objective function. */ t = sqrt(fabs(x[0])); objgrd[0] = sin(t) + 0.5*t*cos(t); t = sqrt(fabs(x[1])); objgrd[1] = sin(t) + 0.5*t*cos(t); } } static void schwefel_confun(Integer *mode, Integer ncnln, Integer n, Integer tdcjsl, const Integer *needc, const double *x, double *c, double *cjsl, Integer nstate, Nag_Comm *comm) { /* Scalars */ double t1, t2; Integer k; if (*mode==0 || *mode==2) { /* Constraint values are required. * Only those for which needc is non-zero need be set. */ for (k = 0; k < ncnln; k++) { if (needc[k]>0) { switch (k) { case 0: c[k] = pow(x[0], 2.0) - pow(x[1], 2.0) + 3.0*x[0]*x[1]; break; case 1: c[k] = cos(pow((x[0]/200.0), 2.0) + (x[1]/100.0)); break; } } } } if (*mode==1 || *mode==2) { /* Constraint derivatives are required. */ for (k = 0; k < ncnln; k++) { switch (k) { case 0: cjsl[0] = 2.0*x[0] + 3.0*x[1]; cjsl[1] = -2.0*x[1] + 3.0*x[0]; break; case 1: t1 = x[0]/200.0; t2 = x[1]/100.0; cjsl[n] = -sin(pow(t1, 2.0) + t2) * (2.0*t1)/200.0; cjsl[n+1] = -sin(pow(t1, 2.0) + t2)/100.0; break; } } } } static void mystart(Integer npts, double quas[], Integer n, Nag_Boolean repeat, const double bl[], const double bu[], Nag_Comm *comm, Integer *mode) { /* Scalars */ Integer i; for (i=0; i