View source for Draft:Gradient-based Repair Differential Evolution

{{AfC submission/draft}}

[[File:Funcionamiendo del algoritmo G-DEmi.gif|thumb|G-DEmi applied to mixed-integer nonlinear programming problems shows its capacity to efficiently explore the discontinuous feasible regions (red lines).]]
The Algorithm G-Demi, an acronym for '''Gradient-based Differential Evolution for Mixed-Integer Nonlinear Programming''', is a variant of the [[differential evolution]] designed to solve mixed-integer nonlinear programming problems  (MINLP).<ref>{{Cite journal|title=Efficiently handling constraints in mixed-integer nonlinear programming problems using gradient-based repair differential evolution|last2=Portilla-Flores|first2=Edgar Alfredo|date=2024-05-31|journal=PeerJ Computer Science|volume=10|pages=e2095|language=en|issn=2376-5992|doi=10.7717/peerj-cs.2095|last3=Mezura-Montes|first3=Efrén|last4=Vega-Alvarado|first4=Eduardo|last5=Calva-Yañez|first5=María Bárbara|last1=Molina-Pérez|first1=Daniel|doi-access=free |pmid=38855217 |pmc=11157599 }}</ref> The presence of both continuous and discrete variables, along with constraints, leads to discontinuous feasible parts of varying sizes in the search space. Traditional evolutionary algorithms face difficulties with these problems due to their insensitivity in handling constraints, leading to the generation of many infeasible solutions. G-DEmi addresses this limitation by integrating a gradient-based repair method within the differential evolution framework. The aim of the repair method is to fix promising infeasible solutions in different subproblems using the gradient information of the constraint set.

== G-DEmi ==
G-DEmi continuously improves a population of candidate solutions through an iterative cycle of generation, evaluation, and selection of trial vectors. In each iteration, new vectors are generated by combining existing solutions. They are evaluated based on their performance and repaired as necessary to satisfy the constraints.

=== Initial Population ===
The initial population <math>\textbf{P}_g</math> is generated by taking random values. For the real variables, random real values are generated, and for the integer variables, random integer values are generated, corresponding to the solution vector <math>[\textbf{x}, \textbf{y}]</math>. Subsequently, the objective function <math>f(\textbf{x}_g)</math> and the degree of constraint violation <math>G(\textbf{x}_g)</math> are evaluated.

=== Mutation and Crossover ===
For each target vector <math>\textbf{x}_{g,i}</math>, a trial vector <math>\textbf{u}_{g,i}</math> is generated using mutation and binomial crossover (<math>rand/1/bin</math>). The integer variables in <math>\textbf{u}_{g,i}</math> are rounded before evaluating the vector in the objective function and constraints.

=== Evaluation and Selection ===
The trial vector is compared with its corresponding target vector, and the better one is selected according to the following feasibility rules:<ref>{{Cite journal|url=https://www.sciencedirect.com/science/article/pii/S0045782599003898|title=An efficient constraint handling method for genetic algorithms|date=2000-06-09|journal=Computer Methods in Applied Mechanics and Engineering|volume=186|issue=2|pages=311–338|access-date=2024-06-12|issn=0045-7825|doi=10.1016/S0045-7825(99)00389-8|last=Deb|first=Kalyanmoy|bibcode=2000CMAME.186..311D }}</ref>

# Between two infeasible solutions, the one with lower constraint violation is preferred.
# If one solution is infeasible and the other one is feasible, the feasible solution is preferred.
# Between two feasible solutions, the one with better objective function value is preferred.

However, If the trial vector fails to improve its target but still has a lower objective function value <math>(f(\textbf{u}_{g,i}) < f(\textbf{x}_{g,i}))</math>, and no other vector of the same subproblem has been repaired in the current generation<math>(\mathbf{y}_{g,i}^\text{u} \notin \mathbf{Y})</math>, this trial vector is repaired.

=== Reparation and Improvement ===
The better solution between the repaired vector and its target vector is passed to the population of the next generation <math>\textbf{P}_{g+1}</math>. Through these steps, G-DEmi generates a new population in each generation.

The following pseudocode illustrates the algorithm:
    '''algorithm''' G-DEmi Framework
    '''input:''' <math>NP</math>, <math>CR</math>, <math>F</math>, <math>k_{max}</math>, <math>T_{min}</math>, <math>Eval_{max}</math>
    '''output:''' The best solution so far
    initialize the population <math>P_g</math>
    evaluate <math>f(x_g)</math> and <math>G(x_g)</math> for each individual in <math>P_g</math>
    <math>Eval \gets 1</math>
    <math>g \gets 1</math>
    '''while''' <math>Eval < Eval_{max}</math> '''do'''
        <math>Y \gets \emptyset</math>
        '''for each''' individual <math>x_{g,i}</math> '''in''' <math>P_g</math> '''do'''
            generate a trial vector <math>u_{g,i}</math>
            round the integer variables in <math>u_{g,i}</math>
            evaluate <math>f(u_{g,i})</math> and <math>G(u_{g,i})</math>
            '''if''' <math>u_{g,i}</math> is better than <math>x_{g,i}</math> '''then'''
                store <math>u_{g,i}</math> into <math>P_{g+1}</math>
            '''elseif''' <math>f(u_{g,i}) < f(x_{g,i})</math> '''and''' <math>y_{g,i}^u \notin Y</math> '''then'''
                repair <math>u_{g,i}</math> 
                evaluate <math>f(u_{g,i})</math> and <math>G(u_{g,i})</math>
                '''if''' <math>u_{g,i}</math> is better than <math>x_{g,i}</math> '''then'''
                    store <math>u_{g,i}</math> into <math>P_{g+1}</math>
                '''end if'''
                <math>Y = Y \cup y_{g,i}^u</math>
            '''end if'''
            update <math>Eval</math>
        '''end for'''
        <math>g \gets g + 1</math>
    '''end while'''

== Gradient-based Repair Method ==
[[File:Reparación en evolución diferencial.png|thumb|Solution Repair in G-DEmi]]
The gradient-based repair method is a crucial component of G-DEmi, designed to address infeasibility in trial vectors generated by differential evolution operators. This method focuses on independently exploring subproblems defined by integer variables. Specifically, to repair a vector with mixed variables <math>[\textbf{x}, \textbf{y}]</math>, only the real variables <math>\textbf{x}</math> are modified while the integer variables <math>\textbf{y}</math> remain fixed.

The method repairs only those trial vectors <math>\textbf{u}</math> that satisfy two conditions: (i) they lost the tournament against their target vectors but have a better objective function value, and (ii) they belong to a subproblem <math>\textbf{y}</math> where no solution has been repaired in the current generation. These two conditions aim to promote the repair of trial vectors with higher potential and ensure that each subproblem is explored independently, avoiding the repair of similar solutions multiple times.

=== Definition of Constraint Violation ===
The constraint violation <math> \mathbf{V}(\mathbf{x}) </math> is defined as a vector that contains the violation degree for each inequality <math> g </math> and equality <math> h </math> constraint in a given problem, for a particular solution vector <math> \mathbf{x} </math>. Parameters <math> n </math> and <math> m </math> denote the number of inequality and equality constraints, respectively, and <math> \epsilon </math> specifies the tolerance for equality constraints. The sign function <math> \text{sgn}(h) </math> preserves the sign of the equality violation.

<math>
\mathbf{V}(\mathbf{x}) = \begin{bmatrix}
\max(g_1(\mathbf{x}), 0) \\
\vdots \\
\max(g_n(\mathbf{x}), 0) \\
\text{sgn}(h_1) \cdot \max(|h_1(\mathbf{x})| - \epsilon, 0) \\
\vdots \\
\text{sgn}(h_m) \cdot \max(|h_m(\mathbf{x})| - \epsilon, 0)
\end{bmatrix}
</math>

=== The Gradient Matrix of Constraints ===
The gradient matrix of these constraints with respect to the <math> N </math>components of <math> \mathbf{x} </math>, denoted as <math> \nabla \mathbf{V}(\mathbf{x}) </math>, is defined as:

<math>
\nabla \mathbf{V}(\mathbf{x}) = \begin{bmatrix}
\frac{\partial g_1(\mathbf{x})}{\partial x_1} & \frac{\partial g_1(\mathbf{x})}{\partial x_2} & \cdots & \frac{\partial g_1(\mathbf{x})}{\partial x_N} \\
\vdots & \vdots & & \vdots \\
\frac{\partial g_n(\mathbf{x})}{\partial x_1} & \frac{\partial g_n(\mathbf{x})}{\partial x_2} & \cdots & \frac{\partial g_n(\mathbf{x})}{\partial x_N} \\
\frac{\partial h_1(\mathbf{x})}{\partial x_1} & \frac{\partial h_1(\mathbf{x})}{\partial x_2} & \cdots & \frac{\partial h_1(\mathbf{x})}{\partial x_N} \\
\vdots & \vdots & & \vdots \\
\frac{\partial h_m(\mathbf{x})}{\partial x_1} & \frac{\partial h_m(\mathbf{x})}{\partial x_2} & \cdots & \frac{\partial h_m(\mathbf{x})}{\partial x_N}
\end{bmatrix}
</math>

=== Finite Difference Approximation ===
The forward finite difference approximation provides an estimate of these derivatives, defined as:

<math>
\nabla \mathbf{V}(\mathbf{x})_{i,j} \approx \frac{f_i(\mathbf{x} + \Delta x \cdot \mathbf{e}_j) - f_i(\mathbf{x})}{\Delta x}
</math>

where <math> \Delta x </math> represents the step size and <math> \mathbf{e}_j </math> is a unitary vector of the same dimension as <math> \mathbf{x} </math>, with a value of 1 for the <math> j </math> component and 0 for the rest.

=== Repair Procedure ===
This repair method aims to transform <math> \mathbf{x} </math> into a feasible solution, which involves adjusting the elements of the vector <math> \mathbf{V}(\mathbf{x}) </math> to zero. Iteratively, a repaired vector <math> \mathbf{x}_{k+1} </math> can be obtained using [[Newton's method|Newton-Raphson's method]] through the following equation, which represents a linear approximation of <math> \mathbf{V}(\mathbf{x}_k) </math> in the direction of the origin:

<math>
\mathbf{x}_{k+1} = \mathbf{x}_k - \nabla \mathbf{V}(\mathbf{x}_k)^{-1} \mathbf{V}(\mathbf{x}_k)
</math>

However, it is common that the number of variables differs from the number of constraints. In this case, the <math> \nabla \mathbf{V}(\mathbf{x}_k) </math> matrix is non-invertible and the [[Moore–Penrose inverse|Moore-Penrose pseudoinverse]] must be used

<math>
\mathbf{x}_{k+1} = \mathbf{x}_k - \nabla \mathbf{V}(\mathbf{x}_k)^+ \mathbf{V}(\mathbf{x}_k)
</math>

Where  <math> \nabla \mathbf{V}(\mathbf{x}_k)^+ </math> represents the pseudoinverse  matrix of the gradient matrix <math> \nabla \mathbf{V}(\mathbf{x}_k) </math>. A computationally efficient way of finding <math> \nabla \mathbf{V}(\mathbf{x}_k)^+ </math> is by employing [[singular value decomposition]].

=== Mixed Variables Repair ===
A mixed trial vector <math> \mathbf{u} = [\mathbf{x}_u^k, \mathbf{y}_u]^T </math> is defined, where only the <math> \mathbf{x}_u^k </math> component is updated during the iterative repair process. As a result, the constraint violation degree vector and the gradient matrix can be defined as:

<math>\mathbf{V}(\mathbf{u}) = \mathbf{V}(\mathbf{x}_u^k, \mathbf{y}_u)</math>

<math> \nabla \mathbf{V}(\mathbf{u}) = \frac{\partial \mathbf{V}(\mathbf{x}_u^k, \mathbf{y}_u)}{\partial \mathbf{x}_u^k} </math>

The repair method follows these steps:
   '''Algorithm:''' Gradient-based repair method
   '''Input:''' <math>\mathbf{u}, k_{max}, T_{min}</math>
   '''Output:''' <math>\mathbf{u}</math>
   Initialize <math>k = 1</math>.
   '''While''' none of the stopping criteria is fulfilled:
     Calculate <math>\mathbf{V}(\mathbf{x}^\text{u}_k, \mathbf{y}^\text{u})</math> 
     Calculate <math>\nabla \mathbf{V}(\mathbf{x}^\text{u}_k, \mathbf{y}^\text{u})</math> 
     Remove zero elements of <math>\mathbf{V}(\mathbf{x}^\text{u}_k, \mathbf{y}^\text{u}) 
      </math> and their corresponding values in <math>\nabla \mathbf{V}(\mathbf{x}^\text{u}_k, 
      \mathbf{y}^\text{u})</math>.
     Calculate the pseudoinverse <math>\nabla \mathbf{V}(\mathbf{x}^\text{u}_k, 
      \mathbf{y}^\text{u})^{+}</math>.
     Calculate <math>\mathbf{x}^\text{u}_{k+1}</math>
     Update <math>\mathbf{x}^\text{u}{k} \leftarrow \mathbf{x}^\text{u}{k+1}</math>.
     Update <math>\mathbf{u} = [\mathbf{x}^\text{u}_k, \mathbf{y}^\text{u}]^\text{T}</math>.
     Increment <math>k \leftarrow k+1</math>.
   '''End While'''
   '''Stopping criteria:'''
   <math>k \geq k_{max}</math>: Maximum number of iterations reached.
   <math>\mathbf{V} = 0</math>: All elements of <math>\mathbf{V}</math> are equal to zero.
   <math>T_u \leq T_{min}</math>: Maximum absolute difference between <math>\mathbf{x}^\text{u} 
      {k+1}</math> and <math>\mathbf{x}^\text{u}{k}</math> is equal to or lower than <math>T_{min} 
      </math>.

== Example of Mixed Variables Repair ==
This repair procedure can be illustrated by the following example. Consider a scenario with one inequality constraint and one equality constraint, as shown below:

<math> g_1(\mathbf{x},y) = {x_1}^2+ {x_2}^2+ {y}^2- 12 \le 0 </math>

<math> h_1(\mathbf{x},y)= x_1+x_2+ {y}- 5.5=0 </math>

Suppose <math> \mathbf{u} = [2 \ 1 \ 1]^T </math> and an equality tolerance <math> \epsilon = 1 \times 10^{-04} </math>. In the first iteration (where <math> k = 1 </math>), <math> \mathbf{x}^\text{u}_1 = [2 \ 1]^T </math> and <math> \mathbf{y}^\text{u} = 1 </math>. Therefore, the vectors <math> \mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u}) </math> and <math> \nabla \mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u}) </math> are computed as follows:<math display="block">
\mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u}) = 
\begin{bmatrix}
\max(-6, 0) \\
-\max(1.5 - \epsilon, 0)
\end{bmatrix}
=
\begin{bmatrix}
0 \\
-1.4999
\end{bmatrix}
</math><math display="block">
\nabla \mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u}) = 
\begin{bmatrix}
\frac{\partial g_1(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u})}{\partial x_1} & \frac{\partial g_1(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u})}{\partial x_2} \\[1ex]
\frac{\partial h_1(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u})}{\partial x_1} & \frac{\partial h_1(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u})}{\partial x_2}
\end{bmatrix}
=
\begin{bmatrix}
4 & 2 \\
1 & 1
\end{bmatrix}
</math>As you can see, only <math> h_1(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})</math> was violated. Therefore, the element of <math> g_1(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})</math> needs to be removed from <math> \mathbf{V}(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})</math> along with its corresponding values in <math> \nabla \mathbf{V}(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})</math>. This leads to <math> \mathbf{V}(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u}) = -1.4999</math>. Then, <math> \nabla \mathbf{V}(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})</math> and its pseudoinverse  <math> \nabla \mathbf{V}(\mathbf{x}^\text{u}_1,\ \mathbf{y}^\text{u})^{+}</math> are computed as follows:<math display="block">
\nabla \mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u}) = 
\begin{bmatrix}
1 & 1
\end{bmatrix}
\Rightarrow
\nabla \mathbf{V}(\mathbf{x}^\text{u}_1, \mathbf{y}^\text{u})^{+} = 
\begin{bmatrix}
0.5 \\
0.5
\end{bmatrix}
</math>Subsequently, the vector <math>\mathbf{x}^\text{u}_2</math> is obtained as follows:<math display="block">
\mathbf{x}^\text{u}_2 = 
\begin{bmatrix}
2 \\
1
\end{bmatrix}
-
\begin{bmatrix}
0.5 \\
0.5
\end{bmatrix}
\begin{bmatrix}
-1.4999
\end{bmatrix}
=
\begin{bmatrix}
2.75 \\
1.75
\end{bmatrix}
</math>The updated vector <math>\mathbf{V}(\mathbf{x}^\text{u}_2, \mathbf{y}^\text{u})</math> results in:<math display="block">
 \mathbf{V}(\mathbf{x}^\text{u}_2, \mathbf{y}^\text{u}) =
 \begin{bmatrix}
\max(-0.3750, 0) \\
\max(0-\epsilon, 0)
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0
\end{bmatrix}
</math>As you can see, the values of <math>(\mathbf{x}^\text{u}_2, \mathbf{y}^\text{u})</math> satisfy all the constraints. Therefore, the trial vector has been successfully repaired, and its new values are <math>\mathbf{u}=[2.75 \ 1.75 \ 1]^\text{T}</math>.

== References ==
<references></references>

== External links ==

<nowiki>
[[Categoría:Combinatorial optimization]]
[[Categoría:Mathematical optimization]]
[[Categoría:Artificial intelligence]]
[[Categoría:Evolutionary computation]]</nowiki>