/* nag_glopt_bnd_pso (e05sac) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #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 monmod(Integer ndim, Integer npar, double x[], const double xb[], double fb, const double xbest[], const double fbest[], const Integer itt[], Nag_Comm *comm, Integer *inform); #ifdef __cplusplus } #endif static void display_result(Integer ndim, const double xb[], const double x_target[], double fb, double f_target, const Integer itt[], Integer inform); static void display_option(const char *optstr, Nag_VariableType optype, Integer ivalue, double rvalue, const char *cvalue); static void get_known_solution(Integer ndim, double x_target[], double *f_target); /* Global constants - set the behaviour of the monitoring function.*/ static const Integer detail_level = 0; static const Integer report_freq = 100; int main(void) { /* This example program demonstrates how to use * nag_glopt_bnd_pso (e05sac) in standard execution, and with a * selection of coupled local minimizers. * The non-default option 'Repeatability = On' is used here, giving * repeatable results. */ /* Scalars */ Integer ndim = 2, npar = 5; Integer exit_status = 0, lcvalue = 17; Integer liopts = 100, lopts = 100; double fb, f_target, rvalue; Integer i, inform, ivalue; /* Arrays */ static double ruser[2] = {-1.0, -1.0}; char cvalue[17], optstr[81]; double *bl = 0, *bu = 0, opts[100], *xb = 0, *x_target = 0; Integer iopts[100], itt[6]; /* Nag Types */ Nag_VariableType optype; NagError fail; Nag_Comm comm; /* Print advisory information.*/ printf("nag_glopt_bnd_pso (e05sac) Example Program Results\n\n"); /* For communication with user-supplied functions: */ comm.user = ruser; printf("Minimization of the Schwefel function.\n\n"); /* Allocate memory.*/ if (!(bl = NAG_ALLOC(ndim, double)) || !(bu = NAG_ALLOC(ndim, double)) || !(xb = NAG_ALLOC(ndim, double)) || !(x_target = NAG_ALLOC(ndim, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Store the known solution of the problem for a comparison.*/ get_known_solution(ndim, x_target, &f_target); /* Set box bounds for problem.*/ for (i = 0; i < ndim; i++) { bl[i] = -500.0; bu[i] = 500.0; } /* Initialize fail structures.*/ INIT_FAIL(fail); /* Initialize the option arrays for nag_glopt_bnd_pso (e05sac) * using nag_glopt_opt_set (e05zkc). */ nag_glopt_opt_set("Initialize = e05sac", iopts, liopts, opts, lopts, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Query some default option values.*/ printf(" Default Option Queries:\n\n"); /* nag_glopt_opt_get (e05zlc). * Option getting routine for nag_glopt_bnd_pso (e05sac). */ nag_glopt_opt_get("Boundary", &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, &fail); if (fail.code == NE_NOERROR) { display_option("Boundary", optype, ivalue, rvalue, cvalue); nag_glopt_opt_get("Maximum Iterations Completed", &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, &fail); } if (fail.code == NE_NOERROR) { display_option("Maximum Iterations Completed", optype, ivalue, rvalue, cvalue); nag_glopt_opt_get("Distance Tolerance", &ivalue, &rvalue, cvalue, lcvalue, &optype, iopts, opts, &fail); } if (fail.code != NE_NOERROR) { printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message); exit_status = 1; goto END; } display_option("Distance Tolerance", optype, ivalue, rvalue, cvalue); /* ------------------------------------------------------------------*/ printf("\n1. Solution without using coupled local minimizer.\n\n"); /* Set various options to non-default values if required.*/ nag_glopt_opt_set("Repeatability = On", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Verify Gradients = Off", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Boundary = Hyperspherical", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Maximum iterations static = 150", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Repulsion Initialize = 30", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Repulsion Finalize = 30", iopts, liopts, opts, lopts, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_glopt_bnd_pso (e05sac). * Global optimization using particle swarm algorithm (PSO), * bound constraints only. */ nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, 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, xb, x_target, fb, f_target, itt, inform); break; case NE_USER_STOP: /* Exit flag set in objfun or monmod and returned in inform.*/ display_result(ndim, xb, x_target, fb, f_target, itt, inform); break; default: /* An error was detected.*/ exit_status = 1; printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message); goto END; } /* ------------------------------------------------------------------*/ printf("2. Solution using coupled local minimizer " "nag_opt_simplex_easy (e04cbc).\n\n"); /* Initialize fail structures.*/ INIT_FAIL(fail); /* Set an objective target.*/ sprintf(optstr, "Target Objective Value = %32.16e", f_target); nag_glopt_opt_set(optstr, iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Target Objective Tolerance = 1.0e-5", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Target Objective Safeguard = 1.0e-8", iopts, liopts, opts, lopts, &fail); /* Set the local minimizer to be nag_opt_simplex_easy (e04cbc) * and set corresponding options. */ if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Minimizer = e04cbc", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Interior Iterations = 10", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Interior Tolerance = 1.0e-4", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Exterior Tolerance = 1.0e-4", iopts, liopts, opts, lopts, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Search for the global optimum.*/ nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, 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, xb, x_target, fb, f_target, itt, inform); break; case NE_USER_STOP: /* Exit flag set in objfun or monmod and returned in inform.*/ display_result(ndim, xb, x_target, fb, f_target, itt, inform); break; default: /* An error was detected.*/ exit_status = 1; printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message); goto END; } /* -----------------------------------------------------------------*/ printf("3. Solution using coupled local minimizer " "nag_opt_conj_grad (e04dgc).\n\n"); /* Initialize fail structures.*/ INIT_FAIL(fail); /* set the local minimizer to be nag_opt_conj_grad (e04dgc) * and set corresponding options. */ nag_glopt_opt_set("Local Minimizer = e04dgc", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Interior Iterations = 5", iopts, liopts, opts, lopts, &fail); if (fail.code == NE_NOERROR) nag_glopt_opt_set("Local Exterior Iterations = 20", iopts, liopts, opts, lopts, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_glopt_opt_set (e05zkc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Search for the global optimum.*/ nag_glopt_bnd_pso(ndim, npar, xb, &fb, bl, bu, objfun_schwefel, 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, xb, x_target, fb, f_target, itt, inform); break; case NE_USER_STOP: /* Exit flag set in objfun or monmod and returned in inform.*/ display_result(ndim, xb, x_target, fb, f_target, itt, inform); break; default: /* An error was detected.*/ exit_status = 1; printf("Error from nag_glopt_bnd_pso (e05sac)\n%s\n", fail.message); goto END; } END: /* Clean up memory.*/ NAG_FREE(bl); NAG_FREE(bu); NAG_FREE(xb); NAG_FREE(x_target); 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; if (comm->user[0] == -1.0) { printf("(User-supplied callback objfun_schwefel, first invocation.)\n"); comm->user[0] = 0.0; } /* 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 monmod(Integer ndim, Integer npar, double x[], const double xb[], double fb, const double xbest[], const double fbest[], const Integer itt[], Nag_Comm *comm, Integer *inform) { Integer i, j; #define X(J, I) x[(J-1)*ndim + (I-1)] #define XBEST(J, I) xbest[(J-1)*ndim + (I-1)] if (comm->user[1] == -1.0) { printf("(User-supplied callback monmod, first invocation.)\n"); comm->user[1] = 0.0; } if (detail_level) { /* Report on the first iteration, and every report_freq iterations.*/ if (itt[0] == 1 || itt[0]%report_freq == 0) { printf("* Locations of particles\n"); for (j = 1; j <= npar; j++) { printf(" * Particle %2ld\n", j); for (i = 1; i <= ndim; i++) printf(" %2ld %13.5f\n", i, X(j, i)); } printf("* Cognitive memory\n"); for (j = 1; j <= npar; j++) { printf(" * Particle %2ld\n", j); printf(" * Best position\n"); for (i = 1; i <= ndim; i++) printf(" %2ld %13.5f\n", i, XBEST(j, i)); printf(" * Function value at best position\n"); printf(" %13.5f\n", fbest[j - 1]); } printf("* Current global optimum candidate\n"); for (i = 1; i <= ndim; i++) printf(" %2ld %13.5f\n", i, xb[i - 1]); printf("* Current global optimum value\n"); printf(" %13.5f\n\n", fb); } } /* If required set *inform<0 to force exit.*/ *inform = 0; #undef XBEST #undef X } 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, const double xb[], const double x_target[], double fb, double f_target, 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\n", "Total particles re-initialized", itt[5]); /* 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; 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 objective optimum : %13.5f\n", f_target); printf(" Achieved objective value : %13.5f\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"); } static void get_known_solution(Integer ndim, double x_target[], double *f_target) { /* Fill in the known solution of multidimensional Schwefel's function. */ Integer i; if (f_target && x_target && ndim > 0) { *f_target = -418.9828872724337*ndim; for (i = 0; i < ndim; i++) x_target[i] = -420.9687463599820; } }