Simulation and Optimizing of 2D mail delivery system with Monte Carlo¶
By Team EpicCool (Zayan, Kabir and Prannaya)
8 and 9 December 2021
Imports¶
In [ ]:
import Pkg;
Pkg.add("Plots")
Pkg.add("Statistics")
Pkg.add("StatsBase")
Pkg.add("WebIO")
Pkg.add("DelimitedFiles")
Pkg.add("PyPlot") #on-the-fly animation using pyplot backend seems less flickers
Pkg.add("Evolutionary")
Pkg.add("Optim")
Pkg.add("NLopt")
Pkg.build("PyCall")
In [1]:
using Plots, Statistics,Random, StatsBase, DelimitedFiles
Rules and policy¶
In this subject, we consider a mail delivery system where mails are delivered on a network of connected stations,
and they can only be delivered from station to station, say via mail service workers on horses.
It starts at $0$ (the source), with two types of mails (A and B. It can be more than two types in a more complicated setup).
Stations $3$(B) and $6$(A) are the respective destinations. After going to the correct destination, the mails are sent to $7$ (the drain).
Here a few basic rules for delivery the mails:
- For each time step, the mails move only one step.
- One station can have at most one mail at each time step.
- The mail is allowed to be delivered to any of the available connected stations.
- The probabilities of going to any station is user-defined
Core functions for delivery algorithm¶
The connectivity matrix describes all connections of the stations to one another.
We indicate 1 in the corresponding row and columns to indicate a connection.
For eg. where 1 connects to 2, the matrix (row,column)->(1,2)=1
At this stage, this is manually input.
The connectivity matrix describing the above node connections that allows for all directions is as follows:
Connectivity matrix for the delivery path (tested)¶
In [2]:
# connectivity matrix
# Vetically, the mail can be delivered back; Horizontally, the mail can only be delivered forward to the destination.
function connect_matrix_2D()
#1st row is source, last row is drain. all entries are connected to source, all exits are connected to drain!
C = zeros(Float64,8,8);
C[1,2] = 1;
C[2,3] = 1;
C[2,5] = 1;
C[3,2] = 1;
C[3,4] = 1;
C[3,6] = 1;
C[4,3] = 1;
C[4,8] = 1;
C[5,2] = 1;
C[5,6] = 1;
C[6,3] = 1;
C[6,5] = 1;
C[6,7] = 1;
C[7,6] = 1;
C[7,8] = 1;
C[8,8] = 1;
return C
end
Out[2]:
In [3]:
function connect_matrix_A()
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = 0.1;
A[2,5] = 0.9;
A[3,2] = 0.1;
A[3,4] = 0.1;
A[3,6] = 0.8;
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = 0.1;
A[5,6] = 0.9;
A[6,3] = 0.1;
A[6,5] = 0.1;
A[6,7] = 0.8;
A[7,8] = 1;
return A
end
Out[3]:
In [4]:
function connect_matrix_B()
B = zeros(Float64,8,8);
B[1,2] = 1;
B[2,3] = 0.9;
B[2,5] = 0.1;
B[3,2] = 0.1;
B[3,4] = 0.8;
B[3,6] = 0.1;
B[4,8] = 1;
B[5,2] = 0.9;
B[5,6] = 0.1;
B[6,3] = 0.8;
B[6,5] = 0.1;
B[6,7] = 0.1;
B[7,6] = 1; #if B mail reaches A(6), mail returns to station 5 with 100% probability.
return B
end
Out[4]:
In [5]:
C = connect_matrix_2D()
Out[5]:
In [6]:
conn_matA = connect_matrix_A()
Out[6]:
In [7]:
conn_matB = connect_matrix_B()
Out[7]:
Function that initializes the matrix of histories (tested)¶
Entries of history matrix, $H$, are initialized to 0.
All mails get ready at position 0.
This matrix stores the trajectory of the mails across time.
Each row corresponds to a mail in the system and each column stores the positions of all the mails at a snapshot in time.
Going across the row tells you the position of a particular mail across time.
In [8]:
function initialize_H(NC,tM)
# tM maximum number of time steps
# NC is the number of mails
return H = zeros(Int16,NC, tM+2); #this tM+2 can be alterd to taste
end
Out[8]:
Function that finds available next positions (tested)¶
Given a mail's position (for example, ip), the following procedures are used to determine the available positions
- Determine the reachable positions by inspecting the corresponding row of the connectivity matrix, $C$
- Determine the occupied positions by inspecting the column of history matrix, $H$, corresponds to the time step.
However, $0$, $L+1$ are excluded.
Recall the position $0$ is the source which hosts the awaiting of mails, and $L+1$ is for the drain of mails, respectively. - Positions available are the positions which are reachable but not occupied. Equivalently, set of reachable subtract set of occupied.
In [9]:
function find_pos_available_k(ic,it,H,C,L,conn_mat)
# ic is the mail number
# it is the time step
# H matrix of histories
# C connectivity matrix
# L is the number of stations
# for a particular mail at a particular time,
aux = size(H);
NC = copy(aux[1]); #NC returns the total number of mails
vec_pos0 = Int64.(zeros(0));
vec_pos = Int64.(zeros(0));
vec_prob0 = Float64.(zeros(0));
vec_prob = Float64.(zeros(0));
ip = copy(H[ic,it]); #obtain current position of the mail from history matrix
for iL = 1:L+2 #check all possible positions, +2 to account for source and drain
if C[ip+1,iL] == 1 #(row,col) check whether this station is connected to station ip (it would be a 1 in the connectivity matrix)
#ip+1 to account for numbering of station, 0th station starts at 1st row.
if conn_mat[ip+1,iL] != 0 #check whether mail type has prob of moving to station iL, if so, store position and corresponding probability
vec_pos0 = append!(vec_pos0, iL-1); #stores all reachable stations
vec_prob0 = append!(vec_prob0, conn_mat[ip+1,iL]); #store prob of respective reachable stations
end
end
end
#conflict resolution
aux = size(vec_pos0) #total number of reachable stations
for ipos = 1:aux[1] #going across all possible stations
s=0;
if (vec_pos0[ipos] != L+1) #if the mail is not outside, which is at L+1
for iNC = 1:NC #going across all mails
if H[iNC,it] == copy(vec_pos0[ipos]); #checks if stations is already occupied by other mails. if so, add 1 to s.
s+=1;
end
end
end
if s==0 #if s!=1, proceed
vec_pos = append!(vec_pos,vec_pos0[ipos]); # if there are no collisions, store particular reachable station in vec_pos which is returned by this function.
vec_prob = append!(vec_prob,vec_prob0[ipos]);
end
end
if !isempty(vec_prob) #if probability vector is not empty, normalise it
vec_prob = vec_prob./sum(vec_prob);#normalise probability here!
end
return vec_pos, vec_prob
end
Out[9]:
In [10]:
H = initialize_H(4,10)
Out[10]:
In [11]:
H[1,1]=7;
H[2,1]=5;
H[3,1]=3;
H
Out[11]:
In [12]:
vec_pos, Avec_prob=find_pos_available_k(2,1,H,C,6,conn_matA)
Out[12]:
Function that delivers a single mail (tested)¶
The mail worker first decides whether to deliver the mail according to the dictated probability. If so, it randomly picks one position from available positions and registers the intention into the corresponding column of history matrix, $H$. Note that, at this stage, the movement is not final due to possible conflicts with other mails.
In [13]:
#Demonstrates how to draw weighted sampling
# Pkg.add("StatsBase") # Only do this once, obviously
# using StatsBase
# items = ["a", 2, 5, "h", "hello", 3]
# weights = [0.1, 0.1, 0.2, 0.2, 0.1, 0.3]
# sample(items, Weights(weights))
# choice_weighted = sample(vec_pos, Weights(Avec_prob))
In [14]:
function move_position_k(ic,it,H,vec_pos, Avec_prob,conn_matA)
#applies probability of moving via prob connectivity matrix A and B
current_pos = copy(H[ic,it]); # current mail position
aux = size(vec_pos);
n_pos = copy(aux[1]);
if n_pos > 0
H[ic,it+1] = sample(vec_pos, Weights(Avec_prob)); #draw randomly from weighted distribution
end
return H
end
Out[14]:
In [15]:
function deliver_mail_k(ic,numA,it,H,C,L,conn_matA,conn_matB)
# numA is number of mails belonging to A. Above numA, all mails belong to B
# if A execute connectivity matrix A
# ic is the mail number
# H matrix of histories
# it time step
# C connectivity matrix
H[ic, it+1] = copy(H[ic, it]);
# vec_pos, vec_prob= find_pos_available_k(ic,it,H,C,L);
if ic <= numA
vec_pos, Avec_prob= find_pos_available_k(ic,it,H,C,L,conn_matA);
H = move_position_k(ic,it,H,vec_pos, Avec_prob,conn_matA)
else
vec_pos, Bvec_prob= find_pos_available_k(ic,it,H,C,L,conn_matB);
H = move_position_k(ic,it,H,vec_pos, Bvec_prob,conn_matB)
end
return H
end
Out[15]:
Check for and resolve conflicts (tested)¶
Given a position in the mail delivery path, following procedure is used to resolve conflicts:
- Determine a list of mails delivering to the same station by inspecting the column of history matrix, $H$, correspond to current time step. For example, given $\text{it}$ (time step) and $\text{iL}$ (position of the station in connected network), the index of entries in $H[:,\text{it}]==\text{iL}$ are the mails which are going into the station iL.
- More than one mails going into the same position constitue a conflict.
- To resolve conflict, we randomly let one mail to move on and make the rest stay put.
- Iterate over all positions in the connected delivery network, rinse and repeat.
In [16]:
function find_resolve_conflicts(H,L,it)
#conflicts only happen for new positions. conflicts from old to new has been resolved in find_pos_available
# it is the time step
# H matrix of histories
# L is the number of positions
Hn = copy(H);
aux = size(H);
NC = aux[1];
# finds list of cars that are in conflict.
for iL = 1:L
list_cars = Int64.(zeros(0));
for ic = 1:NC
if H[ic,it] == iL
list_cars = append!(list_cars, ic);
end
end
aux = size(list_cars);
#in list of cars with conflict, randomly let one car progress and keep the rest stationary.
if aux[1]>1
i_pos = rand(1:aux[1],1,1);
keep_pos = copy(i_pos[1]);
for ipos = 1:aux[1]
Hn[list_cars[ipos], it] = copy(H[list_cars[ipos], it-1]); #the rest are sent back to their prev station
end
Hn[list_cars[keep_pos], it] = copy(H[list_cars[keep_pos], it]); #selected mail proceeds to new position
end
end
return Hn
end
# H = initialize_H(4,4)
# H[1,2] = 1
# H[2,2] = 1
# H[3,2] = 1
# H[4,2] = 1
# H = find_resolve_conflicts(H,L,2)
Out[16]:
Function that delivery all the mails (tested)¶
All mails are to be delivered to the next station, and conflicts are resolved, the movements are final.
In [17]:
function next_step_k(H,it,C,L,numA,conn_matA,conn_matB)
# H matrix of histories
# p_mov probability of moving
# it time step
# C connectivity matrix
aux = size(H);
NC = copy(aux[1]);
for ic = 1:NC
H = deliver_mail_k(ic,numA,it,H,C,L,conn_matA,conn_matB); #
end
H = find_resolve_conflicts(H,L,it+1); #
return H
end
Out[17]:
Example 1¶
In this example, we simulate a delivery process consisting of two kinds of mails (A and B). They are supposed to be sent to two destinations.
In [18]:
NMail = 4; # number of mails
tM = 10; # total time steps
L = 6; # number of sites
numA = 2 # mail A number
numB = NMail - numA # mail B number
C = connect_matrix_2D();
conn_matA = connect_matrix_A()
conn_matB = connect_matrix_B()
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
H
Out[18]:
Measurement¶
Analyze the time spent travelling between station a and station b¶
We determine the time spent for each type of mail (A or B) to get from station a to station b. The distribution is visualized using histogram. Following that, average and standard deviation of time spent are reported. The function itself returns a list of mail numbers corresponding to those which accomplished the journey between station a and station b.|
In [19]:
function time_start_end(H,numA,numB)
# calculate the time taken by each mail from initial station to destination position given the history matrix
# H history matrix
sv = size(H);
max_time = sv[2];
A_first_time_check = zeros(Int16, numA, 1);
B_first_time_check = zeros(Int16, numB, 1);
A_second_time_check = zeros(Int16, numA, 1);
B_second_time_check = zeros(Int16, numB, 1);
init_pos_A = 1;
init_pos_B = 1;
fina_pos_A = 6;
fina_pos_B = 3;
#find first and second time for A mails
for i in 1:numA
aux = findfirst(x -> x == init_pos_A, H[i,:]);
if aux == nothing
aux = 0;
end
A_first_time_check[i] = aux;
aux = findfirst(x -> x == fina_pos_A, H[i,:]);
if aux == nothing
aux = 0;
end
A_second_time_check[i] = aux;
end
#find first and second time for B mails
for i in 1:numB
aux = findfirst(x -> x == init_pos_B, H[numA+i,:]);
if aux == nothing
aux = 0;
end
B_first_time_check[i] = aux;
aux = findfirst(x -> x == fina_pos_B, H[numA+i,:]);
if aux == nothing
aux = 0;
end
B_second_time_check[i] = aux;
end
#compute the difference between first and last time
diff_A = A_second_time_check - A_first_time_check;
diff_B = B_second_time_check - B_first_time_check;
# if the result is 0, it means that it never entered and it should not be counted
# if the results is negative, it means that it never reached the destination, so now we do not count it
time_diff_A = sort(diff_A[findall(x->x>0, diff_A)]);
time_diff_B = sort(diff_B[findall(x->x>0, diff_B)]);
return time_diff_A, time_diff_B
end
Out[19]:
In [20]:
time_diff_A, time_diff_B = time_start_end(H,numA,numB);
display(time_diff_A)
display(time_diff_B)
Evolution Testing¶
Mock Test and Fitness Functions¶
In [21]:
function test_function(x)
NMail = 10; # number of mails
tM = 10; # total time steps
L = 6; # number of sites
numA = 5 # mail A number
numB = NMail - numA # mail B number
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = x[1];
A[2,5] = x[2];
A[3,2] = x[3];
A[3,4] = x[4];
A[3,6] = x[5];
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = x[6];
A[5,6] = x[7];
A[6,3] = x[8];
A[6,5] = x[9];
A[6,7] = x[10];
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
return time_start_end(H,numA,numB)
end
Out[21]:
In [22]:
function fitness_function(x)
NMail = 10; # number of mails
tM = 10; # total time steps
L = 6; # number of sites
numA = 5 # mail A number
numB = NMail - numA # mail B number
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = x[1];
A[2,5] = x[2];
A[3,2] = x[3];
A[3,4] = x[4];
A[3,6] = x[5];
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = x[6];
A[5,6] = x[7];
A[6,3] = x[8];
A[6,5] = x[9];
A[6,7] = x[10];
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
time_diff_A_overall = Int64.(zeros(0));
time_diff_B_overall = Int64.(zeros(0));
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
time_diff_A, time_diff_B = time_start_end(H,numA,numB)
time_diff_A_overall = append!(time_diff_A_overall, time_diff_A)
time_diff_B_overall = append!(time_diff_B_overall, time_diff_B)
end
return sum(time_diff_A_overall) * 1/size(time_diff_A_overall)[1] + sum(time_diff_B_overall) * 1/size(time_diff_B_overall)[1]
end
Out[22]:
Testing Evolutionary (Genetic Algorithms)¶
In [23]:
using Evolutionary
In [24]:
result = Evolutionary.optimize(fitness_function, [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0],
[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],
GA(populationSize = 100, selection = susinv, crossover = discrete, mutation = domainrange(ones(10))),
Evolutionary.Options(iterations=100))
Out[24]:
In [25]:
x = Evolutionary.minimizer(result)
Out[25]:
Node-Based Optimization¶
Utility and Helper Functions + Constants¶
In [26]:
using Optim
In [27]:
function sigmoid(x)
return 1/(1 + exp(-x))
end
Out[27]:
In [28]:
global_x = [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,]
Out[28]:
In [29]:
NMail = 10; # number of mails
tM = 30; # total time steps
L = 6; # number of sites
numA = 5 # mail A number
numB = NMail - numA # mail B number
Out[29]:
Node-Based Fitness Functions¶
In [30]:
#2 values
function fitness_function_node_2(x)
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(x[1]);
A[2,5] = sigmoid(x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[30]:
In [97]:
A = zeros(Float64,8,8);
A[1,2] = 1;
#A[2,3] = sigmoid(x[1]);
#A[2,5] = sigmoid(x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
display(A)
In [31]:
#2 values
function fitness_function_node_3(x)
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(global_x[1]);
A[2,5] = sigmoid(global_x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(x[1]);
A[3,4] = sigmoid(x[2]);
A[3,6] = sigmoid(x[3]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[31]:
In [32]:
function fitness_function_node_5(x)
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(global_x[1]);
A[2,5] = sigmoid(global_x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(x[1]);
A[5,6] = sigmoid(x[2]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[32]:
In [33]:
function fitness_function_node_6(x)
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(global_x[1]);
A[2,5] = sigmoid(global_x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(x[1]);
A[6,5] = sigmoid(x[2]);
A[6,7] = sigmoid(x[3]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[33]:
Testing Functions¶
In [34]:
function test_function_2()
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(global_x[1]);
A[2,5] = sigmoid(global_x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[34]:
In [35]:
function test_function_3(x)
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = x[1];
A[2,5] = x[2];
A[3,2] = x[3];
A[3,4] = x[4];
A[3,6] = x[5];
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = x[6];
A[5,6] = x[7];
A[6,3] = x[8];
A[6,5] = x[9];
A[6,7] = x[10];
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
loss = 0.0
for i = 1:100
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
if(count(x->x==7, H[:,end]) == 0)
loss += 2
else
loss += 1.0/(count(x->x==7, H[:,end]))
end
end
return loss
end
Out[35]:
In [36]:
function individual_test()
A = zeros(Float64,8,8);
A[1,2] = 1;
A[2,3] = sigmoid(global_x[1]);
A[2,5] = sigmoid(global_x[2]);
A[2,3] = A[2,3] * (1.0/(A[2,3] + A[2,5]))
A[2,5] = A[2,5] * (1.0/(A[2,3] + A[2,5]))
A[3,2] = sigmoid(global_x[3]);
A[3,4] = sigmoid(global_x[4]);
A[3,6] = sigmoid(global_x[5]);
A[3,2] = A[3,2] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,4] = A[3,4] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[3,6] = A[3,6] * (1.0/(A[3,2] + A[3,4] + A[3,6]))
A[4,3] = 1; #if A mail reaches B(3), mail returns to station 2 with 100% probability.
A[5,2] = sigmoid(global_x[6]);
A[5,6] = sigmoid(global_x[7]);
A[5,2] = A[5,2] * (1.0/(A[5,2] + A[5,6]))
A[5,6] = A[5,6] * (1.0/(A[5,2] + A[5,6]))
A[6,3] = sigmoid(global_x[8]);
A[6,5] = sigmoid(global_x[9]);
A[6,7] = sigmoid(global_x[10]);
A[6,3] = A[6,3] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,5] = A[6,5] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[6,7] = A[6,7] * (1.0/(A[6,3] + A[6,5] + A[6,7]))
A[7,8] = 1;
C = connect_matrix_2D();
B = connect_matrix_B();
H = initialize_H(NMail,tM);
for it = 1:tM+1
H = next_step_k(H,it,C,L,numA,conn_matA,conn_matB);
end
return H
end
Out[36]:
Optimization Algorithm¶
In [37]:
for i = 1:50
node = rand([2,3,5,6])
if(node == 2)
result = Optim.optimize(fitness_function_node_2, global_x[1:2], SimulatedAnnealing(), Optim.Options(iterations=100))
global_x[1:2] = Optim.minimizer(result)
elseif(node == 3)
result = Optim.optimize(fitness_function_node_3, global_x[3:5], SimulatedAnnealing(), Optim.Options(iterations=100))
global_x[3:5] = Optim.minimizer(result)
elseif(node == 5)
result = Optim.optimize(fitness_function_node_5, global_x[6:7], SimulatedAnnealing(), Optim.Options(iterations=100))
global_x[6:7] = Optim.minimizer(result)
else
result = Optim.optimize(fitness_function_node_6, global_x[8:10], SimulatedAnnealing(), Optim.Options(iterations=100))
global_x[8:10] = Optim.minimizer(result)
end
end
Verification of Results and Comparison with Base Case¶
Predicted Model Loss¶
In [75]:
test_function_2()
Out[75]:
Base Model Loss¶
In [93]:
test_function_3([0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5])
Out[93]:
History Matrix, H of our Predicted Model¶
In [94]:
individual_test()
Out[94]: