_CONSTRAINED OPTIMIZATION_
by Rainer Storn
edited by Bruce Schneier

Listing One

int accept(int dim,float par_vec[],int adapt,int *finish_ptr, float constraint[],int *nocont)
/**C*F************************************************************************
** SRC-FUNCTION   :accept()                                                 **
** LONG_NAME      :accept                                                   **
** AUTHOR         :Dr. Rainer Storn                                         **
** DESCRIPTION    :accept() tests whether a parameter vector par_vec[]      **
**                 falls into the region of acceptability (ROA). If it does **
**                 hit=1 is returned, otherwise hit=0.                      **
** FUNCTIONS      :none                                                     **
** GLOBALS        :none                                                     **
** PARAMETERS     :dim          number of vector elements.                  **
**                 par_vec[]    contains the vector with dim                **
**                              gaussian distributed variables.             **
**                 adapt        control variable.                           **
**                              0 : no adaptation of ROA                    **
**                              1 : adaptation of ROA permitted (shrinkage) **
**                              2 : adaptation of ROA permitted             **
**                                  (relaxation, expansion)                 **
**                 finish_ptr   indicates meeting of requirements.          **
**                              0 : goals not yet reached.                  **
**                              1 : goals has been reached.                 **
**                 constraint[] contains current constraints of ROA.        **
**                 *nocont      contains number of constraints.             **
** PRECONDITIONS  :par_vec[] must contain valid parameter vector.           **
** POSTCONDITIONS :Elements of constraint[] will probably be altered.       **
**                 nocont returns number of constraints to assist           **
**                 printing routine of main().                              **
***C*F*E**********************************************************************/
{
  int   hit, i, n;
  float goal[NARY], z0, z1, x, y;
/*------Set up constraints----------------------------------------------*/
/*------goal[0] : Maximum absolute value of polynomial for x ex [0,1]-----
--------goal[1] : Maximum absolute value of polynomial for x = +/- 1.2--*/
  goal[0] = 1.001;    /* must-constraint */
  goal[1] = 5.9;      /* must-constraint */
  *nocont = 2;        /* two constraints */
/*------Calculate function values and initializations-------------------*/
/*------Passband. Compute maximum magnitude of ordinate value.----------*/
  z0 = 0.;
  for (i=0; i<=100; i++)
  {
    y = 0.0;
    x = -1.0 + (float)i/50;
    for (n=dim-1; n>0; n=n-1)
      {
        y = (y + par_vec[n])*x;
      }
    y = y + par_vec[0];
    if (fabs(y) > z0) z0 = fabs(y);  /* z0 contains maximum magnitude */
  }
/*--------Stopband. Compute ordinate value at the edges x = (+/- 1.2).--*/
/*--------Save the lowest ordinate value.-------------------------------*/
  y  = 0.0;
  x  = 1.2;
  for (n=dim-1; n>0; n=n-1)
  {
    y = (y + par_vec[n])*x;
  }
  y = y + par_vec[0];
  z1 = y;
  y = 0.0;
  x = -1.2;
  for (n=dim-1; n>0; n=n-1)
  {
    y = (y + par_vec[n])*x;
  }
  y = y + par_vec[0];
  if (y < z1) z1 = y;
  hit  = 1;   /* Preset hit-flag to "hit" */
/*------Relax constraints if adapt equals 2.---------------------------------*/
  if (adapt == 2)
  {
    for (i=0; i<=1; i++)  /* Initialize constraints to goal values */
    {
      constraint[i] = goal[i];
    }
/*--Relax must-constraints as required---------------------------------------*/
    if (z0 > constraint[0]) constraint[0] = z0;
    if (z1 < constraint[1]) constraint[1] = z1;
  }
  else if (adapt == 1)    /*--------adapt must-constraints (shrinkage)-------*/
  {
    if (z0 <= constraint[0])
    {
      if (z0 > goal[0]) constraint[0] = z0; /* adapt must-constraint only if */
      else constraint[0] = goal[0];           /* goal is not reached yet.    */
    }
    else
    {
      hit = 0;
    }
    if (z1 >= constraint[1])              /* adapt must-constraint only if */
    {                                     /* goal is not reached yet.      */
      if (z1 < goal[1]) constraint[1] = z1;
      else constraint[1] = goal[1];
    }
    else
    {
      hit = 0;
    }
  }
  else
  {
    if (z0 > constraint[0]) hit = 0;
    if (z1 < constraint[1]) hit = 0;
  }
  *finish_ptr = 0;
  if ((constraint[0] <= goal[0]) && (constraint[1] >= goal[1])) *finish_ptr=1;
  return(hit);
}


Listing Two

200
100
5
0.7
10.00000
10.00000
-6.00000
10.00000
80.00000
3.00100
3.00100
3.00100
3.00100
3.00100


Listing Three

Iteration: 134
---------------
xnominal[0] = 0.993562
xnominal[1] = -0.043474
xnominal[2] = -7.870154
xnominal[3] = 0.050184
xnominal[4] = 7.859358
sigma[0]    = 0.012818
sigma[1]    = 0.009166
sigma[2]    = 0.035221
sigma[3]    = 0.014853
sigma[4]    = 0.040397
constraint[0]    = 1.001000
constraint[1]    = 5.900000
Number of random points = 100
Number of hits          = 3
Yield in percent        = 3.000000