Discrete optimization - Githubissues

Description of problem:

When using the age-structured model, the number of social interactions within the population is controlled using the 'total' interaction matrix $Nc$. This is a 16x16 matrix that, at position X, Y contains the amount of social interaction age bin X has with age bin Y. The total interaction matrix is the sum of the interaction matrix at home, at school, at work, and in other places. In the age-layered model, one can define policies by making linear combinations of the schools, work and 'other' interaction matrices. For instance, a scenario where schools close corresponds to, \begin{equation} N{\text{c,total}} = N{\text{c,home}} + N{\text{c,work}} + N{\text{c,others}}. \end{equation} Opposed to a business-as-usual scenario which corresponds to, \begin{equation} N{\text{c,total}} = N{\text{c,home}} + N{\text{c,schools}} + N{\text{c,work}} + N{\text{c,others}}. \end{equation}

The 'control handle' for our MPC algorithm thus is a discrete scenario in the form of a 16x16 matrix rather than a single 'number'. Naturally, when using an optimization algorithm to find the optimal future policy, a problem arises because the 'sought after' variable is now discrete in nature. One option to find the optimal policy is to simply use brute force, f.i. if the length of the MPC control horizon is 4 and there are for 4 discrete scenarios, then one could nest 4 for loops to calculate all possible combinations of policies and its sum-of-squared errors/economic cost. Of course, this is an unsustainable and hacky way of solving the problem. Ideally, a discrete optimization algorithm is used to avoid unnecessary function evaluations.

Another thing I tried before is assigning each discrete scenario a range of values to which it corresponds. In the code below, 0-1 represents a business-as-usual scenario, 1-2 represents a scenario of strick social distancing along with 48% telework, 70% reduction in 'others' but with the reopening of schools. 2-3 corresponds to the Belgian lockdown. I then used a continuous optimization algorithm (PSO) to find the optimal policies but results were 'moderate'.

# Load interaction matrices
 Nc_home = np.loadtxt("Belgium/BELhome.txt", dtype='f', delimiter='\t')
 Nc_work = np.loadtxt("Belgium/BELwork.txt", dtype='f', delimiter='\t')
 Nc_schools = np.loadtxt("Belgium/BELschools.txt", dtype='f', delimiter='\t')
 Nc_others = np.loadtxt("Belgium/BELothers.txt", dtype='f', delimiter='\t')
 Nc_all = np.loadtxt("Belgium/BELall.txt", dtype='f', delimiter='\t')
 # Use values of thetas to build a list object Ncs containing discrete scenarios
 Ncs=[]
 for i in range(thetas.size):
    if thetas[i]<=1 and thetas[i]>=0:
         Ncs.append(Nc_all)
     elif thetas[i]<=2 and thetas[i]> 1:
         Ncs.append(Nc_home + Nc_schools + 0.01*(1-0.52)*Nc_work + 0.01*(1-0.70)*Nc_others)
     elif thetas[i]<=3 and thetas[i]> 2:
         Ncs.append(Nc_home + 0.01*(1-0.52)*Nc_work + 0.01*(1-0.70)*Nc_others)

Proposed solution: Examine candidate algorithms suited for discrete optimization. An example is GEKKO optimization suite (https://gekko.readthedocs.io/en/latest/). I already tried using it once but this did not work. However, I highly recommend giving GEKKO another try, it looks like a really promising suite.

UGentBiomath / COVID19-Model

Discrete optimization #29