next up previous contents Next: Continuation for one dimensional Up: ALIAS-C++ Previous: Linear algebra

  • La page de J-P. Merlet
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    Subsections



    Optimization

    ALIAS is also able to deal with optimization problem. Let consider a scalar function $F$ of $m$ variables $X=\{x_1,\ldots,x_m\}$, called the optimum function and assume that you are looking for the global minimal or maximal value of $F$ when the unknowns lie in some ranges. Furthermore the unknowns may be submitted to $n_1$ constraints of type $G_i(X)=0$, $n_2$ constraints of type $H_i(X)\le 0$ and $n_3$ constraints of type $K_i(X)\ge 0$. Note that $F, G_i, H_i, K_i$ may have interval coefficients.

    Definition of a minimum and a maximum

    As interval coefficients may appear in the function $F$ we have to define what will be called a minimum or a maximum of $F$. First we assume that there is no interval coefficients in $F$ and denote by $F^m$ the minimal or maximal value of $F$ over a set of ranges defined for $X$ and an accuracy $\epsilon$ with which we want to determine the extremum. The algorithm will return an interval $F_e$ as an approximation of $F^m$ such that for a minimization problem $\underline{F_e} -\epsilon \le F^m \le \underline{F_e}+\epsilon$ and for a maximization problem $\overline{F_e} -\epsilon \le F^m \le \overline{F_e}+\epsilon$. The algorithm will also return a value $X^m$ for $X$ where the extremum occurs. If we deal with a constrained optimization problem we will have:

    where $\epsilon_f$ has a pre-defined value. Note that if we have constraint equation of type $G_i(X)=0$ the result of the optimization problem may be no more guaranteed as the constraint itself $G_i$ may not be verified.

    If there are interval coefficients in the optimum function there is not a unique $F^m$ but according to the value of the coefficient a minimal extremum value $F^m_1$ and a maximal extremum value $F^m_2$. The algorithm will return in the lower bound of $F_e$ an approximation of $F^m_1$ and in the upper bound of $F_e$ an approximation of $F^m_2$ which verify for a minimization problem:

    Note that the width of $F_e$ may now be greater than $\epsilon$.The algorithm will return also two solutions $X^m_1, X^m_2$ for $X$ corresponding respectively to the values of $F^m_1, F^m_2$. If we are dealing with a constrained optimization problem the solutions will verify the above constraint equations.

    If the optimum function has no interval coefficients the algorithm may return no solution if the interval evaluation of the optimum function has a width larger than $\epsilon$. Evidently the algorithm will also return no solution if there is no solution that satisfy all the constraints.

    Methods

    A first method to find all the solutions of an optimum problem is to consider the system of derivative equations and find its root (eventually under the given constraints): this may be done with the solving procedures described in a previous chapter. Hence we have implemented an alternative method which is able to work even if the optimum function has, globally or locally, no derivatives. This method is implemented as a special case of the general solving procedures. The only difference is that the procedure maintains a value for the current optimum: during the bisection process we evaluate the optimum function for each new box and reject those that are not compatible with the current optimum.

    Two types of method enable to solve an optimization problem:

    1. a method which need only a procedure for evaluating $F$,
    2. a method which need a procedure for evaluating $F$ and a procedure for evaluating its gradient,
    Note that with the first method none of the function needs to be differentiable while for the other one not all the functions in the set $F, G, H, K$ must be differentiable: only one function in the set has to be differentiable.

    For the first method we update the current optimum by computing the optimum function value at the middle point of the box. For the method with the gradient a local optimizer based on the steepest descent method is also used if there is only one equation to be minimized or if there are only inequalities constraints. The local optimizer works in 2 steps: first a rough optimization and then a more elaborate procedure if the result of the first step is better than a fixed threshold defined by the global variable ALIAS_Threshold_Optimiser whose default value is 100. Additionally if an extremum $E$ has been already determined the local optimizer (which may computer intensive) is called only if the function value $V$ at the middle point is such that for a maximum $\vert V\vert\ge {\tt ALIAS\_LO\_A} \vert S\vert$ and $\vert V\vert \le {\tt ALIAS\_LO\_B} \vert S\vert$ for a minimum, where ALIAS_LO_A and ALIAS_LO_B are global variables with default value 0.9 and 1.1.


    Implementation

    Preliminary notes:

    Optimization with function evaluation

    The optimization method is implemented as:

     
    int Minimize_Maximize(int m,int n,
          INTEGER_VECTOR &Type_Eq,
          INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), 
          INTERVAL_VECTOR & TheDomain, 
          int Iteration,int Order,
          double epsilon,double epsilonf,double epsilone,
          int Func_Has_Interval,
          INTERVAL Optimum,
          INTERVAL_MATRIX & Solution,
          int (* Simp_Proc)(INTERVAL_VECTOR &));
    
    the arguments being: Thus to minimize a function you have to set its Type_Eq to -2, to maximize it to set its Type_Eq to 2, while if you are looking for both a minimum and a maximum Type_Eq should be set to 10.

    Remember that you may use the 3B method to improve the efficiency of this algorithm (see section 2.3.2) if you have constraint equations.

    In some cases it may be interesting to determine if the minimum and maximum have same sign. This may be done by setting the flag ALIAS_Stop_Sign_Extremum either to:

    With the flag set to 2 the detection that the extremum will have opposite sign is faster.

    Return code

    The procedure will return an integer $k$ that will be identical to the code returned by the procedure
    Solve_General_Interval except for:

    Dealing with inequalities on the same function

    In this version of ALIAS there is no direct way to deal with inequalities that are valid for the same function (e.g. $\alpha \le U(X) \le
\beta$), which mean that you will have to declare two inequalities (which imply that the quantity $U$ will be evaluated twice). But in the previous procedure there is a way to avoid writing two inequalities. In the function evaluation procedure you will just compute $U(X)$ and declare this function as an inequality of type $H(X) \le 0$. After having computed the interval evaluation $U$

    Optimization with function and jacobian evaluation

    The optimization method is implemented as:
     
    int Minimize_Maximize_Gradient(int m,int n,
          INTEGER_VECTOR &Type_Eq,
          INTERVAL_VECTOR (* TheIntervalFunction)(int,int,INTERVAL_VECTOR &), 
          INTERVAL_MATRIX (* Gradient)(int, int,INTERVAL_VECTOR &),
          INTERVAL_VECTOR & TheDomain, 
          int Iteration,int Order,
          double epsilon,double epsilonf,double epsilone,
          int Func_Has_Interval,
          INTERVAL Optimum,
          INTERVAL_MATRIX & Solution,
          int (* Simp_Proc)(INTERVAL_VECTOR &));
    
    the arguments being: Thus to minimize a function you have to set its Type_Eq to -2 and to maximize it to set its Type_Eq to 2.

    Remember that you may use the 3B method to improve the efficiency of this algorithm (see section 2.3.2).

    Note also that a convenient way to write the IntervalFunction and Gradient procedures is to use the possibilities offered by ALIAS-Maple (see the ALIAS-Maple manual).

    Return code

    The procedure will return an integer $k$ that will be identical to the code returned by the procedure
    Solve_General_Gradient_Interval except for:
    Note that this implementation is only a special occurrence of the general solving procedure and thus offer the same possibilities to improve the storage management (see section 2.4).


    Order

    General principle

    During the bisection process new boxes will be created and stored in the list. But we want to order these new boxes so that the procedure will consider first the most promising box. The ordering is based on an evaluation index, the new boxes being stored using an increasing order of the index (the box with the lowest index will be stored first). The flag Order indicate which index is used: indexMAXCONSTRAINTFUNCTION@MAX_CONSTRAINT_FUNCTION


    The variable table

    Assume now that you have chosen a mixed bisection in which the bisection is applied on $n>1$ variables. The procedure will choose the bisected variables using, for example, the smear function. But in some cases it may be interesting to guide the bisection: for example if we know that subsets of the variables have a strong influence on the extremal value of the optimum function it may be interesting to indicate that as soon as the smear function has led to bisecting one variable in a given subset it may be good to bisect also the other variables in the subset. For example consider the following functions:

    \begin{eqnarray*}
&&f_i = u_x +a\cos\theta_i+b\cos(\theta_i+\psi_i) -x_i\\
&&g_i = u_y +a\sin\theta_i+b\sin(\theta_i+\psi_i) -y_i
\end{eqnarray*}

    where $u_x, u_y, a, b, \theta_i, \psi_i$ are unknowns and $x_i, y_i$ are given. Consider now the optimum function $F$:

    \begin{displaymath}
F = \sum_{i =1}^{i =10} f^2_i+g^2_i
\end{displaymath}

    which has 24 unknowns. But clearly each subset $(\theta_i, \psi_i)$ has a strong influence on the minimum of $F$. Hence if one of the $\theta_i$ is bisected it may be interesting to bisect also $\psi_i$. This may be done by setting the flag ALIAS_Guide_Bisection to 1 and using the variables table: for a problem with $m$ unknowns the variables table $V$ is an array of size $m \times m$ and a 1 in $V(i,j)$ indicates that if the variable $i$ is bisected then the variable $j$ should be also bisected. In ALIAS the variables table is implemented under the name ALIAS_Bisecting_Table. It is the responsibility of the user to clear this array and update it as in the following example:
     
    	Resize(ALIAS_Bisecting_Table,24,24);
    	Clear(ALIAS_Bisecting_Table);
    	ALIAS_Bisecting_Table(1,2)=1;
    	ALIAS_Bisecting_Table(2,1)=1;
    	ALIAS_Bisecting_Table(3,4)=1;
    	ALIAS_Bisecting_Table(4,3)=1;
    

    Examples

    Example 1

    Consider the equation: $
\cos(x)+y\cos(y)-0.2=0
$ which defines a curve in the $x-y$ plane on which you want to determine the closest point to the origin when $x,y$ lie in the range $[-\pi,\pi]$. This leads to trying to find the minimum of the function $x^2+y^2$ under the constraint $
\cos(x)+y\cos(y)-0.2=0
$.

    The procedure for the interval evaluation of the 2 functions will be written as:

     
    INTERVAL_VECTOR IntervalTestFunction (int l1,int l2,INTERVAL_VECTOR & in)
    // interval valued test function
    {
    INTERVAL x,y;
    INTERVAL_VECTOR xx(2);
    
    x=in(1);
    y=in(2);
    if(l1==1)
      xx(1)=Cos(x)+y*Sqr(Cos(y))-0.2;
    if(l1<=2 && l2>=2)
      xx(2)=Sqr(x)+Sqr(y);
      return xx;
    }
    
    while the main program may be written as:
     
    int main()
    {
    int Iterations, Dimension,Dimension_Eq,Num,i,j,precision;
    double Accuracy,Accuracy_Variable;
    INTERVAL_MATRIX SolutionList(2,2);
    INTERVAL_VECTOR TestDomain(2),F(2),P(2),H(2);
    INTEGER_VECTOR Type(2);
    INTERVAL Optimum;
    REAL pi;
    
    pi=Constant::Pi;
    
    Dimension_Eq=2;Dimension=2;
    TestDomain(1)=INTERVAL(-pi,pi);TestDomain(2)=INTERVAL(-pi,pi);
    
    cerr << "Number of iteration = "; cin >> Iterations;
    cerr << "Accuracy on Function = "; cin >> Accuracy;
    
    Type(1)=0;Type(2)=-2;
    Accuracy=0;
    
    Num=Minimize_Maximize(Dimension,Dimension_Eq,Type,
           IntervalTestFunction,TestDomain,Iterations,Accuracy_Variable,
           Accuracy,0,Optimum,SolutionList);
    if(Num<0)
      {
        cout << "The procedure has failed, error code:"<<Num<<endl;
        return 0;
      }
    cout<<"Optimum:"<<Optimum<<" obtained at"<<endl;
    for(i=1;i<=Num;i++)
      {
        cout << "x=" << SolutionList(i,1) <<endl;
        cout << "y=" << SolutionList(i,2) <<endl;
      }
      return 0;
    }
    
    The Minimize_Maximize and Minimize_Maximize_Gradient procedures will return the same numerical results but the number of boxes will change. The results obtained for a full bisection,
    the MAX_MIDDLE_FUNCTION_ORDER and according to the accuracy epsilonf and the storage mode (either direct (DSM) or reverse (RSM), see section 2.3.1.2) are presented in the following table (the number of boxes for the Minimize_Maximize_Gradient procedure is indicated in parenthesis):
    epsilonf Minimum $X^m_1$ $X^m_2$ boxes
    0.01, DSM [1.12195667, 1.12195667] [-0.944932,-0.4786] [-0.944932,-0.4786] 76 (36)
    0.01, RSM [1.12195667, 1.12195667] [-0.944932,-0.4786] [-0.944932,-0.4786] 59 (37)
    0.001, DSM [1.1401661,1.1401661] [-0.954903,-0.477835] [-0.954903,-0.477835] 201 (75)
    0.001, RSM [1.1401661,1.1401661] [-0.954903,-0.477835] [-0.954903,-0.477835] 148 (67)
    0.000001, DSM [1.14223267,1.14223267] [-0.957596,-0.474596] [-0.957596,-0.474596] 5031 (2041)
    0.000001, RSM [1.14223267,1.14223267] [-0.957596,-0.474596] [-0.957596,-0.474596] 4590 (2164)
    This example shows clearly the influence of the constraint equation on the determination of the optimum.

    Example 2

    Consider the problem of finding the coordinates $x,y,z$ of a point that lie on the surface $x^3+y^3+z^3+1-(x+y+z+1)^2$, is inside the sphere centered at (-1,-1,0) of radius 1 and is the closest possible to the center of this sphere, with the constraint that $x,y,z$ lie in the range [-2,2]. Thus we have:

    \begin{eqnarray*}
&&F(X)=(x+1)^2+(y+1)^2+z^2\\
&&G(X)=x^3+y^3+z^3+1-(x+y+z+1)^2=0\\
&&H(X)=(x+1)^2+(y+1)^2+z^2-1 \le 0
\end{eqnarray*}

    With epsilonf=0.0001 we find out that the point is located at at (-0.747,-0.747,0.086059) which is well inside the sphere and that the minimal distance is 0.13529.

    We may also compute the minimal distance not to a point but to a line segment, for example defined by $x \in [0.9,1.1]$, $y=-1$, $z=0$. In that case the optimum function in the evaluation procedure may be defined as:

     
    Sqr(x+INTERVAL(0.9,1.1))+Sqr(y+1)+Sqr(z)
    
    and with epsilonf=0.0001 the algorithm will return that the minimal distance lie in the range [0.0907,0.1925].


    next up previous contents Next: Continuation for one dimensional Up: ALIAS-C++ Previous: Linear algebra
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    jean-pierre merlet
    2018-07-25