On the 17th of April from 18:00 to 20:00, We will be upgrading GitLab. During this time the service will be unavailable, apologies for any inconvenience
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> during the in-class session some of the confusion was caused by code issues. Comments relevant to the code and notebooks as-used in the Friday in-class session are provided in boxes like this.</p></div>
%% Cell type:markdown id: tags:
## Introduction
_Note: part of the background material for this project was already available in [Chapter 5.11 of the textbook](https://mude.citg.tudelft.nl/2023/book/optimization/project.html)._
* We showed in the previous notebook how to use MILP to solve the Road Network Design (RND) Problem
* You saw that due to a large number of binary variables it takes long to reach a low gap
* Think about larger problems, like the road network of Amsterdam or Shanghai, and it will be even harder!
* Here we show how a metaheuristic such a the genetic algorithm can be used to find good (not necessarily optimal) solutions for the problem in potentially less time
%% Cell type:markdown id: tags:
## Libraries
To run this notebook you need to have installed the following packages:
Pandas
Numpy
Matplotlib
Gurobipy
PyMOO
Luckily, they're all part of this weeks environment!
%% Cell type:code id: tags:
``` python
importpandasaspd
importnumpyasnp
importgurobipyasgp
importmatplotlib.pyplotasplt
#if you did not install pymoo yet, run the following in your Anaconda Prompt: conda install -c anaconda autograd
# Genetic algorithm dependencies. We are importing the pymoo functions that are imporant for applying GA (the package can also apply other methods)
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> the intention with the commented lines above was to illustrate possible crossover methods you could use, but unfortunately the pymoo documentation was not clear enough to indicate how they could have been used. The methods <code>SinglePointCrossover</code> and <code>TwoPointCrossover</code> (n<em>not</em> listed above!) could have been easily used, as these methods require no additional keyword arguments and would have worked “out of the box;” they were also illustrated with examples in the online textbook. We apologize for this oversight on our part.</p></div>
As we discussed, it is challenging to use MILP for large-scale NDPs. Therefore, in this assignment, we’re going to use a genetic algorithm to address this problem.
Genetic Algorithms (GAs) are powerful optimization techniques inspired by the process of natural evolution. They have gained prominence in solving complex problems across various fields, ranging from engineering and economics to artificial intelligence. Here, we give a brief overview of Genetic Algorithms, highlighting their fundamental principles, components, and applications in solving optimization problems.
At the heart of a Genetic Algorithm are populations of potential solutions, represented as individuals or chromosomes. These individuals evolve over generations to improve their fitness with respect to the given optimization objective.
Basic Components of a Genetic Algorithm:
***Population**: A collection of individuals representing potential solutions to the problem.
***Objective Function** (or fitness function): Quantifies the quality of each individual’s solution with respect to the optimization objective.
***Selection**: Individuals are chosen based on their fitness to serve as parents for the next generation.
***Crossover**: Genetic material from parents is combined to create offspring with potentially improved traits.
***Mutation**: Random changes are introduced in offspring to maintain diversity and explore new solution spaces.
***Replacement**: New offspring replace some of the least fit individuals in the population.
***Termination Criteria**: Conditions under which the algorithm stops, e.g., a maximum number of generations or satisfactory fitness level.
%% Cell type:markdown id: tags:
### GA steps
Reminding you about the GA steps …
* Initialization (start): A population of random individuals is generated to start the algorithm.
* Evaluation (fitness assessment): The fitness function assesses the quality of each individual’s solution.
* Selection: Individuals with higher fitness have a higher chance of being selected as parents.
* Crossover: Genetic material (part of the solutions) from selected parents is combined to create offspring.
* Mutation: Random changes are introduced to some offspring to maintain diversity.
* Replacement: New offspring replace some individuals in the population.
* Termination (end): The algorithm stops when a termination criterion is met.
%% Cell type:markdown id: tags:
### PyMOO
PyMOO is a Python library that provides a comprehensive and easy-to-use framework for multi-objective optimization (MOO). For this case, we are going to deal with only one objective; nevertheless, this is an useful tool if you have more objectives. In addition, PyMOO easily allows us to define our optimization problem by specifying the objectives, constraints, and decision variables.
%% Cell type:markdown id: tags:
## Problem definition and formulation
Here is the problem formulation as presented in the previous Jupyter notebook
& \sum_{j \in N; (i,j) \in A}{ x_{ijs}} - \sum_{j \in N; (j,i) \in A}{ x_{jis}} = d_{is} \quad \forall i \in N, \forall s \in D \\
& y_{ij} \in \{0, 1\}\quad \forall (i,j) \in A \\
& x_{ij} \geq 0 \quad \forall (i,j) \in A \\
& x_{ijs} \geq 0 \quad \forall (i,j) \in A, \forall s \in D \\
\end{align}$$
You can check the NDP notebook for details. But here we deal with it slightly differently to be able to use GA. Essentially, we break down the problem into two sub-problems: 1) the traffic assignment (TA) problem: the route choices of the drivers, and the 2) the road network design problem (NDP): where we select which links should be upgraded. We solve the problem by iteratively going between the Traffic assignment and the Design Problem. The idea is for the GA to move to better networks as generations pass which are evaluated by the traffic assignment process that you have learned.
We use Gurobi to solve the Traffic Assignment sub-problems, which provide us with the objective function (or fitness function within the context of GA) value of the decision problem (which will be dealt with using GA). This is usually referred to as the iterative-optimization-assignment method since we iteratively improve the objective function value of the NDP using the assignment problem.
So let's see how that works.
%% Cell type:markdown id: tags:
### The network design sub-problem
The network desing is where we use the genetic algorithm. As explained before, GA uses a population of solutions and iteratively improves this population to evolve to new generations of populations with a better objective function value (being that minimization or maximization). In this problem, the decision variables are links for capacity expansion and the objective function value is the total system travel time that we want to minimize.
Where the values of $x_{ij}$ are not decision variables anymore, they will be obtained from solving the Traffic Assignment problem with Gurobi which evaluates the travel times on the network. This part of the problem will not be solved mathematically anymore, the $y_{ij}$ variables are decided by the genetic algorithm through the process you learned.
%% Cell type:markdown id: tags:
### The traffic assignment sub-problem
This is just part of the original NDP that assigns traffic to the network based on a set of given capacity values, which are defined based on the values of the DP decision variables (links selected for capacity expansion). The main difference (and the advantage) here is that by separating the binary decision variables, instead of a mixed integer programming problem, which are hard to solve, here we have a quadratic programming problem with continuous decision variables, which will be transformed to a linear problem that Gurobi can solve very fast.
& \sum_{j \in N; (i,j) \in A}{ x_{ijs}} - \sum_{j \in N; (j,i) \in A}{ x_{jis}} = d_{is} \quad \forall i \in N, \forall s \in D \\
& x_{ij} \geq 0 \quad \forall (i,j) \in A \\
& x_{ijs} \geq 0 \quad \forall (i,j) \in A, \forall s \in D \\
\end{align}
Where the values of $y_{ij}$ are constant and are defined by GA.
%% Cell type:markdown id: tags:
### Summarizing
The following is a diagram that shows what you are finally doing to solve the same problem but with a meta-heuristic approach:
%% Cell type:markdown id: tags:
## Data preprocessing
Our data preprocessing steps are similar to the previous notebook. We use some networks from the well-known transportation networks for benchmarking repository as well as a small toy network for case studies of NDPs. the following functions read data from this repository and perform data preprocessing to have the input and the parameters required for our case studies.
# create node-destination demand matrix (not a regular OD!)
demand={(n,d):0forninnodesfordindestinations}
forrinorigins:
forsindestinations:
if (r,s)inods_data:
demand[r,s]=ods_data[r,s]
forsindestinations:
demand[s,s]=-sum(demand[j,s]forjinorigins)
returndemand
```
%% Cell type:markdown id: tags:
Now that we have the required functions for reading and processing the data, let's define some problem parameters and prepare the input.
%% Cell type:markdown id: tags:
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> the variables <code>extension_max_no</code> and <code>timelimit</code> were both defined here, but never used. <code>extension_max_no</code> is similar to the variable in notebook A with the same name; here it is superceded below by <code>Budget</code>.</p></div>
%% Cell type:code id: tags:
``` python
# define parameters, case study (network) list and the directory where their files are
extension_factor=1.5# capacity after extension (1.5 means add 50%)
extension_max_no=20# the number of links to add capacity to (simplified way of reprsenting a budget for investing)
#it's the same to say that it's exactly this number of that this number is the max, that's because every investment brings travel time benefits
#even if just one car circulates.
timelimit=300# seconds
beta=2# parameter to use in link travel time function (explained later)
G,pos=network_visualization(link_flow=fftts,coordinates_path=coordinates_path)# the network we create here will be used later for further visualizations!
```
%% Cell type:markdown id: tags:
Now we are ready to build our models!
%% Cell type:markdown id: tags:
### Modeling and solving the traffic assignment sub-problem with Gurobi
In this section we build a Gurobi model to solve the Traffic Assignment sub-problems. The decision variables, objective function, and the constraints of this problem were described before.
Here we wrap the code in a function so that we can use it later within the GA.
# In the GA part the only variables are the binary decision variables, don't forget that the traffic assignment that
# produces the travel time on the network is done in the evaluation of the solution
g=sum(decision_vars)-self.budget
out["F"]=total_travel_time
out["G"]=g
```
%% Cell type:markdown id: tags:
Now, let's initiate an instance of the problem based on the problem class we defined, and initiate the GA with its parameters. Note that depending on the problem size and the number of feasible links, you might need larger values for population and generation size to achieve good results or even feasible results. Of course this increases the computation times.
%% Cell type:markdown id: tags:
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> population size <code>pop_size</code> was 10 originally. If you change this, you will see different results. This is problem-dependent!</p></div>
%% Cell type:markdown id: tags:
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>This:</b><code>Budget</code> is the way the number of links was selected (unlike notebook A, which used <code>extension_max_no</code>). </p></div>
%% Cell type:markdown id: tags:
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>This:</b><code>Budget</code> is the way the number of links was selected (unlike notebook A, which used <code>extension_max_no</code>). The initial value of 76 was trivial, and the solution converged quickly. Numbers between 10 and 40 would have produced <em>much</em> more interesting results (see solution explanation).</p></div>
%% Cell type:code id: tags:
``` python
Budget=76
pop_size=10
# initiate an instance of the problem with max number of selected links as budget constraint
problem=NDP(budget=Budget)
# initiate the GA with parameters appropriate for binary variables
method=GA(pop_size=pop_size,
sampling=BinaryRandomSampling(),
mutation=BitflipMutation(),
crossover=HalfUniformCrossover()
)
```
%% Cell type:markdown id: tags:
Now we are ready to minimize the NDP problem using the GA method we defined.
%% Cell type:markdown id: tags:
<divstyle="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> termination is set here as a keyword argument (see note above).</p></div>
%% Cell type:code id: tags:
``` python
opt_results=minimize(problem,
method,
termination=("time","00:05:00"),#5 minute maximum computation time
seed=7,
save_history=True,
verbose=True,
)
print("Best Objective Function value: %s"%opt_results.F)
print("Constraint violation: %s"%opt_results.CV)
print("Best solution found: %s"%opt_results.X)
#To better interpret the results, this is the legend:
#n_gen: Number of generations
#n_eval: Number of function evaluations
#cv_min: Minimum constraint violation
#cv_avg: Average constraint violation
#f_avg: Average objective function value
#f_min: Minimum objective function value
```
%% Cell type:markdown id: tags:
### Convergence curve
Let's first define some functions (to use later) to get the results and plot them.
%% Cell type:code id: tags:
``` python
defget_results(opt_results):
number_of_individuals=[]# The number of individuals in each generation
optimal_values_along_generations=[]# The optimal value found in each generation
forgeneration_statusinopt_results.history:
# retrieve the optimum from the algorithm
optimum=generation_status.opt
# filter out only the feasible solutions and append and objective space values
