/* nag_glopt_nlp_pso (e05sbc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim, const double x[], double *objf, double vecout[], Integer nstate, Nag_Comm *comm); static void NAG_CALL confun(Integer *mode, Integer ncon, Integer ndim, Integer tdcj, const Integer needc[], const double x[], double c[], double cjac[], Integer nstate, Nag_Comm *comm); static void NAG_CALL monmod(Integer ndim, Integer ncon, Integer npar, double x[], const double xb[], double fb, const double cb[], const double xbest[], const double fbest[], const double cbest[], const Integer itt[], Nag_Comm *comm, Integer *inform); static void display_option(const char *optstr, Nag_VariableType optype, Integer ivalue, double rvalue, const char *cvalue); static void display_result(Integer ndim, Integer ncon, const double xb[], double fb, const double cb[], const Integer itt[], Integer inform); #ifdef __cplusplus } #endif /* Global constants.*/ /* Set the behaviour of the monitoring function.*/ static const Integer detail_level = 0; static const Integer report_freq = 100; /* Known solution for a comparison.*/ static const double f_target = -731.70709230672696; static const double c_scale[] = { 2490.0, 750000.0, 0.1 }; static const double c_target[] = { 0.0, 0.0, 0.0 }; static const double x_target[] = { -394.1470221120988, -433.48214189947606 }; int main(void) { /* This example program demonstrates how to use * nag_glopt_nlp_pso (e05sbc) in standard execution, and with * nag_opt_nlp (e04ucc) as a coupled local minimizer. * The non-default option 'Repeatability = On' is used here, giving * repeatable results. */ /* Scalars */ Integer ncon = 3, ndim = 2, npar = 20; Integer exit_status = 0, lcvalue = 17; Integer liopts = 100, lopts = 100; double fb, rvalue; Integer i, inform, ivalue; /* Arrays */ char cvalue[17], optstr[81]; double opts[100], *bl=0, *bu=0, *cb=0, *xb=0; double *cbest=0, *fbest=0, *xbest=0; Integer iopts[100], itt[7]; /* NAG types */ Nag_VariableType optype; Nag_Comm *comm; NagError fail; /* Print advisory information.*/ printf("nag_glopt_nlp_pso (e05sbc) Example Program Results\n\n"); printf("Minimization of the Schwefel function.\n"); printf("Subject to one linear and two nonlinear constraints.\n\n"); /* Allocate memory for arrays.*/ if (!(bl = NAG_ALLOC(ndim+ncon, double)) || !(bu = NAG_ALLOC(ndim+ncon, double)) || !(cb = NAG_ALLOC(ncon, double)) || !(cbest = NAG_ALLOC(ncon*npar, double)) || !(fbest = NAG_ALLOC(npar, double)) || !(xb = NAG_ALLOC(ndim, double)) || !(xbest = NAG_ALLOC(ndim*npar, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i=0; i < npar; i++) fbest[i] = 0.0; for (i=0; i < npar*ndim; i++) xbest[i] = 0.0; for (i=0; i < npar*ncon; i++) cbest[i] = 0.0; /* Set problem specific values.*/ /* Set box bounds.*/ for (i = 0; i < ndim; i++) { bl[i] = -500.0; bu[i] = 500.0; } /* Set constraint bounds.*/ bl[ndim] = -1.0e6; bl[ndim+1] = -1.0; bl[ndim+2] = -0.9; bu[ndim] = 10.0; bu[ndim+1] = 5.0e5; bu[ndim+2] = 0.9; /* Initialize NagError structure.*/ INIT_FAIL(fail); /* Initialize comm.*/ comm = NAGCOMM_NULL; /* Initialize the option arrays for nag_glopt_nlp_pso (e05sbc) * using nag_glopt_opt_set (e05zkc). */ nag_glopt_opt_set("Initialize = e05sbc", iopts, liopts, opts, lopts, NAGERR_DEFAULT); /* Query some default option values.*/ printf(" Default Option Queries:\n\n"); /* nag_glopt_opt_get (e05zlc). * Option getting routine for nag_glopt_nlp_pso (e05sbc). */ ivalue = 0; rvalue = 0.0; optype = 0; strcpy(optstr, "Constraint Norm"); nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, NAGERR_DEFAULT); display_option(optstr, optype, ivalue, rvalue, cvalue); strcpy(optstr, "Maximum Iterations Completed"); nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, NAGERR_DEFAULT); display_option(optstr, optype, ivalue, rvalue, cvalue); strcpy(optstr, "Distance Tolerance"); nag_glopt_opt_get(optstr, &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, NAGERR_DEFAULT); display_option(optstr, optype, ivalue, rvalue, cvalue); /* ------------------------------------------------------------------*/ printf("\n1. Solution without using coupled local minimizer.\n\n"); /* Set various options to non-default values if required.*/ sprintf(optstr, "Distance Tolerance = %32.16e", rvalue * 0.1); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Constraint Tolerance = 1.0e-4", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Constraint Norm = Euclidean", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Repeatability = On", iopts, liopts, opts, lopts, NAGERR_DEFAULT); sprintf(optstr, "Target Objective Value = %32.16e", f_target); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Target Objective Tolerance = 1.0e-4", iopts, liopts, opts, lopts, NAGERR_DEFAULT); /* nag_glopt_nlp_pso (e05sbc). Global optimization using particle swarm algorithm (PSO), comprehensive. */ nag_glopt_nlp_pso(ndim, ncon, npar, xb, &fb, cb, bl, bu, xbest, fbest, cbest, objfun_schwefel, confun, monmod, iopts, opts, comm, itt, &inform, &fail); /* It is essential to test fail.code on exit.*/ switch (fail.code) { case NE_NOERROR: case NW_FAST_SOLUTION: case NW_SOLUTION_NOT_GUARANTEED: /* No errors, best found solution at xb returned in fb.*/ display_result(ndim, ncon, xb, fb, cb, itt, inform); break; case NE_USER_STOP: /* Exit flag set in objfun, confun or monmod and returned in inform.*/ display_result(ndim, ncon, xb, fb, cb, itt, inform); break; default: /* An error was detected.*/ exit_status = 1; printf("Error from nag_glopt_nlp_pso (e05sbc)\n%s\n", fail.message); goto END; } /* ------------------------------------------------------------------*/ printf("2. Solution using coupled local minimizer nag_opt_nlp (e04ucc).\n\n"); /* Set the local minimizer to be nag_opt_nlp (e04ucc) and set corresponding * options. */ nag_glopt_opt_set("Local Minimizer = e04ucc", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Local Interior Major Iterations = 15", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Local Interior Minor Iterations = 5", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Local Exterior Major Iterations = 50", iopts, liopts, opts, lopts, NAGERR_DEFAULT); nag_glopt_opt_set("Local Exterior Minor Iterations = 15", iopts, liopts, opts, lopts, NAGERR_DEFAULT); /* Query the option Distance Tolerance*/ nag_glopt_opt_get("Distance Tolerance", &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, NAGERR_DEFAULT); /* Adjust Distance Tolerance dependent upon its current value*/ sprintf(optstr, "Distance Tolerance = %32.16e", rvalue*10.0); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT); sprintf(optstr, "Local Interior Tolerance = %32.16e", rvalue); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT); sprintf(optstr, "Local Exterior Tolerance = %32.16e", rvalue * 1.0e-4); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, NAGERR_DEFAULT); /* nag_glopt_nlp_pso (e05sbc). Global optimization using particle swarm algorithm (PSO), comprehensive. */ nag_glopt_nlp_pso(ndim, ncon, npar, xb, &fb, cb, bl, bu, xbest, fbest, cbest, objfun_schwefel, confun, monmod, iopts, opts, comm, itt, &inform, &fail); /* It is essential to test fail.code on exit.*/ switch (fail.code) { case NE_NOERROR: case NW_FAST_SOLUTION: case NW_SOLUTION_NOT_GUARANTEED: case NW_NOT_FEASIBLE: /* nag_glopt_nlp_pso (e05sbc) encountered no errors during * operation, and will have returned the best solution found. */ display_result(ndim, ncon, xb, fb, cb, itt, inform); break; case NE_USER_STOP: /* Exit flag set in objfun, confun or monmod and returned in inform.*/ display_result(ndim, ncon, xb, fb, cb, itt, inform); break; default: /* An error was detected.*/ exit_status = 1; printf("Error from nag_glopt_nlp_pso (e05sbc)\n%s\n", fail.message); goto END; } END: /* Clean up memory.*/ if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); if (cb) NAG_FREE(cb); if (cbest) NAG_FREE(cbest); if (fbest) NAG_FREE(fbest); if (xb) NAG_FREE(xb); if (xbest) NAG_FREE(xbest); return exit_status; } static void NAG_CALL objfun_schwefel(Integer *mode, Integer ndim, const double x[], double *objf, double vecout[], Integer nstate, Nag_Comm *comm) { /* Objective function routine returning the schwefel function and * its gradient. */ Nag_Boolean evalobjf, evalobjg; Integer i; /* Test nstate to indicate what stage of computation has been reached.*/ switch (nstate) { case 2: /* objfun is called for the very first time. */ break; case 1: /* objfun is called on entry to a NAG local minimizer. */ break; default: /* This will be the normal value of nstate. */ ; } /* Test mode to determine whether to calculate objf and/or objgrd.*/ evalobjf = Nag_FALSE; evalobjg = Nag_FALSE; switch (*mode) { case 0: case 5: /* Only the value of the objective function is needed.*/ evalobjf = Nag_TRUE; break; case 1: case 6: /* Only the values of the ndim gradients are required.*/ evalobjg = Nag_TRUE; break; case 2: case 7: /* Both the objective function and the ndim gradients are required.*/ evalobjf = Nag_TRUE; evalobjg = Nag_TRUE; } if (evalobjf) { /* Evaluate the objective function.*/ *objf = 0.0; for (i = 0; i < ndim; i++) *objf += x[i]*sin(sqrt(fabs(x[i]))); } if (evalobjg) { /* Calculate the gradient of the objective function, */ /* and return the result in vecout.*/ for (i = 0; i < ndim; i++) { vecout[i] = sqrt(fabs(x[i])); vecout[i] = sin(vecout[i]) + 0.5*vecout[i]*cos(vecout[i]); } } } static void NAG_CALL confun(Integer *mode, Integer ncon, Integer ndim, Integer tdcj, const Integer needc[], const double x[], double c[], double cjac[], Integer nstate, Nag_Comm *comm) { /* Supplies constraints. * cjac[(k-1)*tdcj + (i-1)] corresponds to dc[k]/dx[i] * for k=1,...,ncon and i=1,...,ndim. */ Integer k; Nag_Boolean evalc, evalcjac; /* Test nstate to determine whether the local minimizer is being called * for the first time from a new start point. */ if (nstate == 1) { /* Set any constant elements of the Jacobian matrix.*/ cjac[0] = 3.0; cjac[1] = -2.0; } /* mode: are constraints, derivatives, or both are required? */ evalc = (*mode == 0 || *mode == 2) ? Nag_TRUE : Nag_FALSE; evalcjac = (*mode == 1 || *mode == 2) ? Nag_TRUE : Nag_FALSE; for (k = 0; k < ncon; k++) { if(needc[k] <= 0) continue; if (evalc == Nag_TRUE) { /* Constraint values are required. * Only those for which needc is non-zero need be set. */ switch (k) { case 0: c[k] = 3.0*x[0] - 2.0*x[1]; break; case 1: c[k] = x[0]*x[0] - x[1]*x[1] + 3.0*x[0]*x[1]; break; case 2: c[k] = cos(pow((x[0]/200.0), 2) + (x[1]/100.0)); break; default: c[k] = 0.0; } } if (evalcjac == Nag_TRUE) { /* Constraint derivatives (cjac) are required.*/ switch (k) { case 0: /* Constant derivatives set when nstate=1 remain throughout * the local minimization. */ break; case 1: /* If the constraint derivatives are known and are readily * calculated, populate cjac when required. */ cjac[k*tdcj] = 2.0*x[0] + 3.0*x[1]; cjac[k*tdcj+1] = -2.0*x[1] + 3.0*x[0]; break; default: /* Any elements of cjac left unaltered will be approximated * using finite differences when required. */ ; } } } /* If an immediate exit is required, return *mode<0 */ *mode = 0; } static void NAG_CALL monmod(Integer ndim, Integer ncon, Integer npar, double x[], const double xb[], double fb, const double cb[], const double xbest[], const double fbest[], const double cbest[], const Integer itt[], Nag_Comm *comm, Integer *inform) { Integer i; if (detail_level) { /* Report on the first iteration, and every report_freq iterations.*/ if (itt[0] == 1 || itt[0]%report_freq == 0) { printf("* Current global optimum candidate\n"); for (i=0; i < ndim; i++) printf(" %2ld %12.5f\n", i+1, xb[i]); printf("* Current global optimum value\n"); printf(" %13.5f\n", fb); printf("* Current constraint violations\n"); for (i = 0; i < ncon; i++) printf(" %2ld %12.5f\n", i+1, cb[i]); printf("\n"); } } /* If required set *inform<0 to force exit.*/ *inform = 0; } static void display_option(const char *optstr, Nag_VariableType optype, Integer ivalue, double rvalue, const char *cvalue) { /* Subroutine to query optype and print the appropriate option values.*/ switch (optype) { case Nag_Integer: printf("%-38s: %13ld\n", optstr, ivalue); break; case Nag_Real: printf("%-38s: %13.4f\n", optstr, rvalue); break; case Nag_Character: printf("%-38s: %13s\n", optstr, cvalue); break; case Nag_Integer_Additional: printf("%-38s: %13ld %16s\n", optstr, ivalue, cvalue); break; case Nag_Real_Additional: printf("%-38s: %13.4f %16s\n", optstr, rvalue, cvalue); break; default:; } } static void display_result(Integer ndim, Integer ncon, const double xb[], double fb, const double cb[], const Integer itt[], Integer inform) { /* Display final results in comparison to known global optimum.*/ Integer i; /* Display final counters.*/ printf(" Algorithm Statistics\n"); printf(" --------------------\n"); printf("%-38s: %13ld\n", "Total complete iterations", itt[0]); printf("%-38s: %13ld\n", "Complete iterations since improvement", itt[1]); printf("%-38s: %13ld\n", "Total particles converged to xb", itt[2]); printf("%-38s: %13ld\n", "Total improvements to global optimum", itt[3]); printf("%-38s: %13ld\n", "Total function evaluations", itt[4]); printf("%-38s: %13ld\n", "Total particles re-initialized", itt[5]); printf("%-38s: %13ld\n\n", "Total constraints violated", itt[6]); /* Display why finalization occurred.*/ switch (inform) { case 1: printf("Solution Status : Target value achieved\n"); break; case 2: printf("Solution Status : Minimum swarm standard deviation obtained\n"); break; case 3: printf("Solution Status : Sufficient number of particles converged\n"); break; case 4: printf("Solution Status : Maximum static iterations attained\n"); break; case 5: printf("Solution Status : Maximum complete iterations attained\n"); break; case 6: printf("Solution Status : Maximum function evaluations exceeded\n"); break; case 7: printf("Solution Status : Feasible point located\n"); break; default: if (inform<0) printf("Solution Status : User termination, inform = %16ld\n", inform); else printf("Solution Status : Termination, an error has been detected\n"); break; } /* Display final objective value and location.*/ printf(" Known constrained objective optimum : %13.3f\n", f_target); printf(" Achieved objective value : %13.3f\n\n", fb); printf(" Comparison between the known optimum and the achieved solution.\n"); printf(" x_target xb\n"); for (i=0; i < ndim; i++) printf(" %2ld %12.2f %12.2f\n", i+1, x_target[i], xb[i]); printf("\n"); if (ncon > 0) { printf(" Comparison between scaled constraint violations.\n"); printf(" c_target cb\n"); for (i = 0; i < ncon; i++) printf(" %2ld %12.5f %12.5f\n", i+1, c_target[i]/c_scale[i], cb[i]/c_scale[i]); printf("\n"); } }