11.2. Tilburg area development case

Python code for the Tilburg area development case (Chapter 9.1). Python code for the shopping mall design problem revisited non-linear example (Chapter 7.2).

Note

Press the rocket symbol on the top right of the page to make the page interactive and play with the code!

Import Packages

This script is fairly similar to the non-linear shopping mall example. Only the preference functions for objective 1 and 2 are changed, together with the weights.

Note that the non-linear preference curves are created by using an interpolation function called pchip_interpolation.

import micropip
await micropip.install("urllib3 ")
await micropip.install("requests")

import sys
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import pchip_interpolate

from genetic_algorithm_pfm import GeneticAlgorithm

Set Weights for Different Objectives

Set weights for the different objectives.

w1 = 1/6
w2 = 1/6
w3 = 1/6
w4 = 1/6
w5 = 1/6
w6 = 1/6

def objective_p1(h1, h2, h3, o, c, s):
    """
    Objective to maximize the profit preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([80, 130, 180], [0, 30, 100], (h1))


def objective_p2(h1, h2, h3, o, c, s):
    """
    Objective to maximize the amount of affordable houses preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([60, 120, 140], [0, 30, 100], (h2))

def objective_p3(h1, h2, h3, o, c, s):
    """
    Objective to maximize the profit preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([60, 70, 80], [0, 50, 100], (h3))


def objective_p4(h1, h2, h3, o, c, s):
    """
    Objective to maximize the amount of affordable houses preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([28, 29, 40], [100, 10, 0], (o))


def objective_p5(h1, h2, h3, o, c, s):
    """
    Objective to maximize the profit preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([30, 68, 80], [100, 30, 0], (c))


def objective_p6(h1, h2, h3, o, c, s):
    """
    Objective to maximize the amount of affordable houses preference.

    :param x1: 1st design variable
    :param x2: 2nd design variable
    :param x3: 3rd design variable
    :param x4: 4th design variable
    :param x5: 5th design variable
    :param x6: 6th design variable
    """
    return pchip_interpolate([30, 35, 45], [0, 70, 100], (s))

def objective(variables):
    """
    Objective function that is fed to the GA. Calles the separate preference functions that are declared above.

    :param variables: array with design variable values per member of the population. Can be split by using array
    slicing
    :return: 1D-array with aggregated preference scores for the members of the population.
    """
    # extract 1D design variable arrays from full 'variables' array
    h1 = variables[:, 0]
    h2 = variables[:, 1]
    h3 = variables[:, 2]
    o  = variables[:, 3]
    c  = variables[:, 4]
    s  = variables[:, 5]

    # calculate the preference scores
    p_1 = objective_p1(h1, h2, h3, o, c, s)
    p_2 = objective_p2(h1, h2, h3, o, c, s)
    p_3 = objective_p3(h1, h2, h3, o, c, s)
    p_4 = objective_p4(h1, h2, h3, o, c, s)
    p_5 = objective_p5(h1, h2, h3, o, c, s)
    p_6 = objective_p6(h1, h2, h3, o, c, s)

    # aggregate preference scores and return this to the GA
    return [w1, w2, w3, w4, w5, w6], [p_1, p_2, p_3, p_4, p_5, p_6]

Define Constraints and Bounds

Before we can run the optimization, we finally need to define the constraints and bounds.

def constraint_1(variables):
    """Constraint that checks if the sum of the areas x1 and x2 is not higher than 10,000 m2.

    :param variables: ndarray of n-by-m, with n the population size of the GA and m the number of variables.
    :return: list with scores of the constraint
    """
    h1 = variables[:, 0]
    h2 = variables[:, 1]
    h3 = variables[:, 2]
    o  = variables[:, 3] * 1000
    c  = variables[:, 4] * 1000
    s  = variables[:, 5] / 100

    return .1 * (100 * h1) + .17 * (80 * h2) + .2 * (120 * h3) + .1 * o + .25 * (40 * h1 + 35 * h2 + 28 * h3) + .33 * c + 100500 * s + 40 * h1 + 40 * h2 + 28 * h3 + .3 * (o + (40 * h1 + 35 * h2 + 28 * h3) + c ) - 100500  # < 0


def constraint_2(variables):
    """Constraint that checks if the sum of the areas x1 and x2 is not lower than 3,000 m2.

    :param variables: ndarray of n-by-m, with n the population size of the GA and m the number of variables.
    :return: list with scores of the constraint
    """
    h1 = variables[:, 0]
    h2 = variables[:, 1]
    h3 = variables[:, 2]
    o  = variables[:, 3] * 1000
    c  = variables[:, 4] * 1000
    s  = variables[:, 5] / 100

    return 85 * h1 + 85 * h2 + 75 * h3 - 100500 * s  # < 0

# define list with constraints
cons = [['ineq', constraint_1]]

# set bounds for all variables
b1 = [80, 180]  # x1
b2 = [60, 140]  # x2
b3 = [60, 80]  # x3
b4 = [28, 40]  # x4
b5 = [30, 80]  # x5
b6 = [30, 45]  # x6
bounds = [b1, b2, b3, b4, b5, b6]

Graphical Output

Setup the graphical output.

# create arrays for plotting continuous preference curves
c1 = np.linspace(80, 180)
c2 = np.linspace(60, 140)
c3 = np.linspace(60, 80)
c4 = np.linspace(28, 40)
c5 = np.linspace(30, 80)
c6 = np.linspace(30, 45)

# calculate the preference functions
p1 = pchip_interpolate([80, 155, 180], [0, 30, 100], (c1))
p2 = pchip_interpolate([60, 120, 140], [0, 30, 100], (c2))
p3 = pchip_interpolate([60, 70, 80], [0, 50, 100], (c3))
p4 = pchip_interpolate([28, 29, 40], [100, 10, 0], (c4))
p5 = pchip_interpolate([30, 68, 80], [0, 30, 100], (c5))
p6 = pchip_interpolate([30, 35, 45], [0, 70, 100], (c6))

# create figure that plots all preference curves and the preference scores of the returned results of the GA
fig = plt.figure(figsize=((10,15)))

font1 = {'size':20}
font2 = {'size':15}

plt.rcParams['font.size'] = '12'
plt.rcParams['savefig.dpi'] = 300

ax1 = fig.add_subplot(3, 2, 1)
ax1.plot(c1, p1, label='Preference curve', color='black')
ax1.set_xlim((80, 180))
ax1.set_ylim((0, 102))
ax1.set_title('Municipality')
ax1.set_xlabel('Number of student houses [-]')
ax1.set_ylabel('Preference score')
ax1.grid()
ax1.legend()
ax1.grid(linestyle = '--')

ax2 = fig.add_subplot(3, 2, 2)
ax2.plot(c2, p2, label='Preference curve', color='black')
ax2.set_xlim((60, 140))
ax2.set_ylim((0, 102))
ax2.set_title('Municipality')
ax2.set_xlabel('Number of starter houses [-]')
ax2.set_ylabel('Preference score')
ax2.grid()
ax2.legend()
ax2.grid(linestyle = '--')

ax3 = fig.add_subplot(3, 2, 3)
ax3.plot(c3, p3, label='Preference curve', color='black')
ax3.set_xlim((60, 80))
ax3.set_ylim((0, 102))
ax3.set_title('Municipality')
ax3.set_xlabel('Number of remaining houses [-]')
ax3.set_ylabel('Preference score')
ax3.grid()
ax3.legend()
ax3.grid(linestyle = '--')

#fig = plt.figure()
ax4 = fig.add_subplot(3, 2, 4)
ax4.plot(c4, p4, label='Preference curve', color='black')
ax4.set_xlim((28, 40))
ax4.set_ylim((0, 102))
ax4.set_title('Municipality')
ax4.set_xlabel('Amount of offices [m2x1000]')
ax4.set_ylabel('Preference score')
ax4.grid()
ax4.legend()
ax4.grid(linestyle = '--')


ax5 = fig.add_subplot(3, 2, 5)
ax5.plot(c5, p5, label='Preference curve', color='black')
ax5.set_xlim((30, 80))
ax5.set_ylim((0, 102))
ax5.set_title('Municipality')
ax5.set_xlabel('Amount of cultural facilities [m2x1000]')
ax5.set_ylabel('Preference score')
ax5.grid()
ax5.legend()
ax5.grid(linestyle = '--')

#fig = plt.figure()
ax6 = fig.add_subplot(3, 2, 6)
ax6.plot(c6, p6, label='Preference curve', color='black')
ax6.set_xlim((30, 45))
ax6.set_ylim((0, 102))
ax6.set_title('Municipality')
ax6.set_xlabel('Open space ratio [-]')
ax6.set_ylabel('Preference score')
ax6.grid()
ax6.legend()
ax6.grid(linestyle = '--')

ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()
ax5.legend()
ax6.legend()

fig.tight_layout()

Optimization

Now we have everything for the optimization, we can run it. For more information about the different options to configure the GA, see the docstring of GeneticAlgorithm (via help()) or chapter 4 of the reader.

# We run the optimization with two paradigms
paradigm = ['minmax', 'a_fine']
marker = ['o', '*']
colours = ['orange', 'green']

# Define the figure and axes before the loop
fig = plt.figure(figsize=(12, 8))

# Creating six subplots for the six preference scores
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 3)
ax4 = fig.add_subplot(2, 3, 4)
ax5 = fig.add_subplot(2, 3, 5)
ax6 = fig.add_subplot(2, 3, 6)

for i in range(2):
    # Dictionary with parameter settings for the GA run with the IMAP solver
    options = {
        'n_bits': 8,
        'n_iter': 400,
        'n_pop': 500,
        'r_cross': 0.8,
        'max_stall': 8,
        'aggregation': paradigm[i],  # minmax or a_fine
        'var_type': 'real'
    }

    # Run the GA and print its result
    print(f'Run GA with {paradigm[i]}')
    ga = GeneticAlgorithm(objective=objective, constraints=cons, bounds=bounds, options=options)
    score_IMAP, design_variables_IMAP, _ = ga.run()

    # Print the optimal result in a readable format
    print(f'Optimal result for h1 = {round(design_variables_IMAP[0], 2)} houses and ' 
          f'h2 = {round(design_variables_IMAP[1], 2)} houses and ' 
          f'h3 = {round(design_variables_IMAP[2], 2)} houses and ' 
          f'o = {round(design_variables_IMAP[3], 2)} m² offices and ' 
          f'c = {round(design_variables_IMAP[4], 2)} m² cultural facilities and ' 
          f's = {round(design_variables_IMAP[5], 2)} open space ratio')

    # Calculate individual preference scores for the results
    c1_res = design_variables_IMAP[0]
    p1_res = pchip_interpolate([80, 155, 180], [0, 30, 100], c1_res)

    c2_res = design_variables_IMAP[1]
    p2_res = pchip_interpolate([60, 120, 140], [0, 30, 100], c2_res)

    c3_res = design_variables_IMAP[2]
    p3_res = pchip_interpolate([60, 70, 80], [0, 50, 100], c3_res)

    c4_res = design_variables_IMAP[3]
    p4_res = pchip_interpolate([28, 29, 40], [100, 10, 0], c4_res)

    c5_res = design_variables_IMAP[4]
    p5_res = pchip_interpolate([30, 68, 80], [0, 30, 100], c5_res)

    c6_res = design_variables_IMAP[5]
    p6_res = pchip_interpolate([30, 35, 45], [0, 70, 100], c6_res)

    # Debugging prints to check calculated values
    print(f"c1_res: {c1_res}, p1_res: {p1_res}")
    print(f"c2_res: {c2_res}, p2_res: {p2_res}")
    print(f"c3_res: {c3_res}, p3_res: {p3_res}")
    print(f"c4_res: {c4_res}, p4_res: {p4_res}")
    print(f"c5_res: {c5_res}, p5_res: {p5_res}")
    print(f"c6_res: {c6_res}, p6_res: {p6_res}")

    # Plot the preference curves 
    if i == 0:  # Only plot curves once
        # create arrays for plotting continuous preference curves
        c1 = np.linspace(80, 180)
        c2 = np.linspace(60, 140)
        c3 = np.linspace(60, 80)
        c4 = np.linspace(28, 40)
        c5 = np.linspace(30, 80)
        c6 = np.linspace(30, 45)

        # calculate the preference functions
        p1 = pchip_interpolate([80, 155, 180], [0, 30, 100], (c1))
        p2 = pchip_interpolate([60, 120, 140], [0, 30, 100], (c2))
        p3 = pchip_interpolate([60, 70, 80], [0, 50, 100], (c3))
        p4 = pchip_interpolate([28, 29, 40], [100, 10, 0], (c4))
        p5 = pchip_interpolate([30, 68, 80], [0, 30, 100], (c5))
        p6 = pchip_interpolate([30, 35, 45], [0, 70, 100], (c6))

        # Plot each preference curve on the respective subplot (assuming ax1 to ax6 are defined)
        ax1.plot(c1, p1, label='Preference curve', color='black')
        ax1.set_xlim((80, 180))
        ax1.set_ylim((0, 102))
        ax1.set_title('Houses type 1')
        ax1.set_xlabel('Number of houses [-]')
        ax1.set_ylabel('Preference score')
        ax1.grid(linestyle='--')

        ax2.plot(c2, p2, label='Preference curve', color='black')
        ax2.set_xlim((60, 140))
        ax2.set_ylim((0, 102))
        ax2.set_title('Houses type 2')
        ax2.set_xlabel('Number of houses [-]')
        ax2.set_ylabel('Preference score')
        ax2.grid(linestyle='--')

        ax3.plot(c3, p3, label='Preference curve', color='black')
        ax3.set_xlim((60, 80))
        ax3.set_ylim((0, 102))
        ax3.set_title('Houses type 3')
        ax3.set_xlabel('Number of houses [-]')
        ax3.set_ylabel('Preference score')
        ax3.grid(linestyle='--')

        ax4.plot(c4, p4, label='Preference curve', color='black')
        ax4.set_xlim((28, 40))
        ax4.set_ylim((0, 102))
        ax4.set_title('Office space')
        ax4.set_xlabel('m2 of office space')
        ax4.set_ylabel('Preference score')
        ax4.grid(linestyle='--')

        ax5.plot(c5, p5, label='Preference curve', color='black')
        ax5.set_xlim((30, 80))
        ax5.set_ylim((0, 102))
        ax5.set_title('Cultural facilities')
        ax5.set_xlabel('m2 of cultural facilities')
        ax5.set_ylabel('Preference score')
        ax5.grid(linestyle='--')

        ax6.plot(c6, p6, label='Preference curve', color='black')
        ax6.set_xlim((30, 45))
        ax6.set_ylim((0, 102))
        ax6.set_title('Open space ratio')
        ax6.set_xlabel('Open space ratio [-]')
        ax6.set_ylabel('Preference score')
        ax6.grid(linestyle='--')

    # Plot the results on the preference curve subplots
    ax1.scatter(c1_res, p1_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])
    ax2.scatter(c2_res, p2_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])
    ax3.scatter(c3_res, p3_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])
    ax4.scatter(c4_res, p4_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])
    ax5.scatter(c5_res, p5_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])
    ax6.scatter(c6_res, p6_res, label='Optimal solution ' + paradigm[i], color=colours[i], marker=marker[i])

# Add legends and adjust layout after plotting all data points
ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()
ax5.legend()
ax6.legend()

# Set titles and labels for the subplots
ax1.set_title('Houses (h1)')
ax1.set_xlabel('Number of houses')
ax1.set_ylabel('Preference score')

ax2.set_title('Houses (h2)')
ax2.set_xlabel('Number of houses')
ax2.set_ylabel('Preference score')

ax3.set_title('Houses (h3)')
ax3.set_xlabel('Number of houses')
ax3.set_ylabel('Preference score')

ax4.set_title('Offices (o)')
ax4.set_xlabel('m² of offices')
ax4.set_ylabel('Preference score')

ax5.set_title('Cultural facilities (c)')
ax5.set_xlabel('m² of cultural facilities')
ax5.set_ylabel('Preference score')

ax6.set_title('Open space ratio (s)')
ax6.set_xlabel('Open space ratio')
ax6.set_ylabel('Preference score')

# Adjust the layout
fig.tight_layout()

# Display the plot
plt.show()

Run GA with minmax
The type of aggregation is set to minmax
Generation   Best score   Mean             Max stall    Diversity    Number of non-feasible results
0            12.7913      734000004.3408   0            0.003        367         
1            11.6111      152000012.1091   0            0.062        76          
2            10.3516      78000011.2888    0            0.123        39          
3            9.2537       26000010.6148    0            0.15         13          
4            5.0175       18000009.723     0            0.157        9           
5            3.4554       22000006.3103    0            0.16         11          
6            3.1002       42000004.4105    0            0.155        21          
7            3.1002       28000003.9718    1            0.16         14          
8            3.0138       22000003.8417    0            0.159        11          
9            3.0138       32000003.8256    1            0.157        16          
10           3.0138       22000003.5345    2            0.159        11          
11           3.0138       28000003.6841    3            0.16         14          
12           3.0138       24000003.6805    4            0.16         12          
13           3.0138       30000003.5307    5            0.155        15          
14           3.0138       38000003.533     6            0.157        19          
15           3.0138       30000003.5249    7            0.159        15          
16           3.0138       24000003.8584    8            0.157        12          
Stopped at gen 16
Execution time was 0.7476 seconds
Optimal result for h1 = 174.14 houses and h2 = 136.56 houses and h3 = 77.58 houses and o = 28.09 m² offices and c = 35.08 m² cultural facilities and s = 37.44 open space ratio
c1_res: 174.140625, p1_res: 79.43010330200194
c2_res: 136.5625, p2_res: 84.97751235961914
c3_res: 77.578125, p3_res: 87.890625
c4_res: 28.09375, p4_res: 90.29370614871809
c5_res: 35.078125, p5_res: 0.5220550998300448
c6_res: 37.44140625, p6_res: 81.91731280004475
Run GA with a_fine
The type of aggregation is set to a_fine
Generation   Best score   Mean             Max stall    Diversity    Number of non-feasible results
No initial starting point for the optimization with the a-fine-aggregator is given. A random population is generated.
0            -100.0       703999988.7618   1            0.004        352         
1            -86.1859     45999962.1966    1            0.062        23          
2            -100.0       51999933.4575    1            0.14         26          
3            -100.0       55999922.9766    1            0.147        28          
4            -100.0       29999921.4447    1            0.149        15          
5            -100.0       27999924.1158    1            0.153        14          
6            -96.7572     43999918.2421    1            0.153        22          
7            -94.9529     49999917.2364    1            0.156        25          
8            -96.9928     43999913.6289    1            0.154        22          
9            -100.0       45999910.4171    1            0.154        23          
10           -96.8457     51999912.5267    1            0.157        26          
11           -96.2197     37999911.2798    1            0.156        19          
12           -96.4365     59999913.4184    2            0.156        30          
13           -93.9386     75999917.3539    3            0.157        38          
14           -100.0       55999908.6982    4            0.158        28          
15           -94.3125     45999913.6609    5            0.158        23          
16           -93.227      53999916.2058    6            0.16         27          
17           -94.1088     55999915.8756    7            0.158        28          
18           -93.1786     57999915.6331    8            0.157        29          
Stopped at gen 18
Execution time was 1.0976 seconds
Optimal result for h1 = 179.61 houses and h2 = 139.69 houses and h3 = 79.92 houses and o = 28.0 m² offices and c = 30.0 m² cultural facilities and s = 41.19 open space ratio
c1_res: 179.609375, p1_res: 98.66712570190427
c2_res: 139.6875, p2_res: 98.6671257019043
c3_res: 79.921875, p3_res: 99.609375
c4_res: 28.0, p4_res: 100.0
c5_res: 30.0, p5_res: 0.0
c6_res: 41.19140625, p6_res: 95.041226417246

Results

Now we have the results, that can be seen in the above figures. First, the resulting design variables are plotted into the solution space. Secondly, we can plot the preference functions together with the results of the optimizations.