This is used in the paper Flexible cable routing framework for wind farm collection system optimization.

Preamble

import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

import numpy as np
import polars as pl

from optiwindnet.db import modelv2
from optiwindnet.db.storagev2 import G_from_routeset
from optiwindnet.svg import svgplot

plt.style.use('jupyter_dark')

%config InlineBackend.figure_formats = ['svg']
plt.rcParams['svg.fonttype'] = 'none'

Setup

db = modelv2.open_database('data/paperdb.v2.sqlite')

Selecting actual routesets

Ns_q = db.NodeSet.select(
    lambda ns:
    ns.name[0] != '!'
)
Ns_q.count()

69

entries = []
for ns in Ns_q:
    for rs in ns.RouteSets:
        entries.append((
            ns.name,
            rs.id,
            ns.T,
            rs.capacity,
            rs.detextra,
            rs.length,
            rs.creator,
            rs.misc.get('gap', None),
        ))
len(entries)

1464

Make DataFrame

schema_actual = dict(
    name=str,
    id=int,
    T=int,
    capacity=int,
    detextra=float,
    length=float,
    creator=str,
    gap=float,
)

df_actual = pl.DataFrame(entries, schema=schema_actual, orient='row')

df_actual

shape: (1_464, 8)

name	id	T	capacity	detextra	length	creator	gap
str	i64	i64	i64	f64	f64	str	f64
"Gode Wind 1"	1	55	5	0.000847	63626.868227	"MILP.pyomo.cplex"	0.004949
"Gode Wind 1"	201	55	2	0.002243	114710.370228	"MILP.pyomo.cplex"	-2.2204e-16
"Gode Wind 1"	5	55	8	null	53268.425754	"MILP.pyomo.cplex"	0.004997
"Gode Wind 1"	4155	55	8	null	54675.524766	"baselines.hgs"	null
"Gode Wind 1"	3617	55	7	null	56703.616107	"baselines.hgs"	null
…	…	…	…	…	…	…	…
"Cazzaro-2022"	21629	50	9	0.0	8.146533	"MILP.pyomo.cplex"	0.005
"Cazzaro-2022"	21613	50	3	-0.007245	14.37826	"baselines.hgs"	null
"Cazzaro-2022"	21623	50	2	-0.018452	19.059407	"MILP.pyomo.cplex"	0.003381
"Cazzaro-2022"	21630	50	10	0.0	7.876699	"MILP.pyomo.cplex"	0.00496
"Cazzaro-2022"	21617	50	7	0.0	8.784742	"baselines.hgs"	null

df_actual['name'].value_counts()

shape: (69, 2)

name	count
str	u32
"Saint-Nazaire"	22
"Saint-Brieuc"	22
"Seagreen"	22
"Rudong H8"	22
"DanTysk"	22
…	…
"Nysted"	22
"Walney 1"	22
"Sandbank"	22
"Anholt"	22
"Belwind"	8

df_actual['gap'].max()

0.019817266284630097

df_actual['detextra'].max()

0.07397519014347398

group ⟨heuri, exact⟩ pairs

pairs_actual = []
for ns in Ns_q:
    dfn = df_actual.filter(pl.col('name') == ns.name)
    for k in range(dfn['capacity'].min(), dfn['capacity'].max() + 1):
        dfk = dfn.filter(pl.col('capacity') == k)
        if dfk.height >= 2:
            mask = pl.col('creator') == 'baselines.hgs'
            heuri = dfk.filter(mask).bottom_k(1, by='length')[0]
            exact = dfk.filter(~mask).bottom_k(1, by='gap')[0]
            if exact['gap'].item() <= 0.02:
                pairs_actual.append((db.RouteSet[heuri['id'].item()],
                                     db.RouteSet[exact['id'].item()]))
len(pairs_actual)

732

diffs = []
for rs_hgs, rs_bnc in pairs_actual:
    diffs.append(rs_hgs.length/rs_bnc.length - 1)

for (heuri, exact), diff in zip(pairs_actual, diffs):
    if diff >= 0.07:
        print(heuri.id, exact.id, exact.nodes.name)

21311 22196 Sandbank
21545 21925 Moray West.1_OSS

small_diffs = []
for (heuri, exact), diff in zip(pairs_actual, diffs):
    if diff <= 0.01:
        small_diffs.append((heuri.id, exact.id))
len(small_diffs)

541

rs = db.RouteSet[21311]
print(rs.detextra)
svgplot(G_from_routeset(rs))

0.07397519014347398

rs = db.RouteSet[22196]
print(rs.detextra)
svgplot(G_from_routeset(rs))

0.0013995437575948788

rs = db.RouteSet[21545]
print(rs.detextra)
svgplot(G_from_routeset(rs))

0.02038470064164777

rs = db.RouteSet[21925]
print(rs.detextra)
svgplot(G_from_routeset(rs))

None

plt.hist(diffs, bins=21);

HGS optimizer

Select HGS-CVRP entries

method_hgs = db.Method.select(
    lambda m:
    m.funname == 'hgs_cvrp'
    and m.options['timeLimit'] == 120
)
method_hgs.count()

2

for m in method_hgs:
    print(m.funname, m.options)

hgs_cvrp {'complete': False, 'lambda_': 40, 'mu': 25, 'nbClose': 5, 'nbElite': 4, 'nbGranular': 20, 'nbIter': 20000, 'scale': 10000.0, 'seed': 0, 'targetFeasible': 0.2, 'timeLimit': 120.0, 'useSwapStar': True}
hgs_cvrp {'complete': False, 'lambda_': 40, 'mu': 25, 'nbClose': 5, 'nbElite': 4, 'nbGranular': 20, 'nbIter': 20000, 'seed': 20241113, 'targetFeasible': 0.2, 'timeLimit': 120, 'useSwapStar': True}

synthetic instances (a single capacity per NodeSet)

pairs_synth = []

# hgs instances that were first solved by cplex
hgs_sub2 = db.RouteSet.select(
    lambda rs:
    rs.method in method_hgs
    and 'solver_details' in rs.misc
    and 'incumbent_id' in rs.misc['solver_details']
    and rs.valid is None
    # and rs.nodes.name != 'Ormonde'
    and rs.nodes.name[0] == '!'
    and rs.T >= 50
)
# keep only instances for which the exact solution is within 2% gap
for rs_hgs in hgs_sub2:
    rs_bnc = db.RouteSet[rs_hgs.misc['solver_details']['incumbent_id']]
    # if rs_bnc.misc['bound']/rs_bnc.length >= 0.98:
    if rs_bnc.misc['gap'] <= 0.02:
        pairs_synth.append((rs_hgs, rs_bnc))
hgs_sub2.count(), len(pairs_synth)

(2346, 2346)

# hgs instances that were solved by cplex afterward
hgs_sub2_ = db.RouteSet.select(
    lambda rs:
    rs.method in method_hgs
    and 'solver_details' not in rs.misc
    and rs.valid is None
    and rs.nodes.name[0] == '!'
    and rs.T >= 50
)
# keep only instances for which the exact solution is within 2% gap
for rs_hgs in hgs_sub2_:
    if len(rs_hgs.nodes.RouteSets) == 1:
        # no MILP counterpart
        continue
    cplex = gurobi = None
    for rs2 in rs_hgs.nodes.RouteSets:
        if rs2 is rs_hgs:
            continue
        if rs2.creator == 'MILP.pyomo.cplex':
            cplex = rs2
        elif rs2.creator == 'MILP.pyomo.gurobi':
            gurobi = rs2
    if gurobi and gurobi.misc['gap'] <= 0.02:
        pairs_synth.append((rs_hgs, gurobi))
    elif cplex and cplex.misc['gap'] <= 0.02:
        pairs_synth.append((rs_hgs, cplex))
hgs_sub2_.count(), len(pairs_synth)

(4033, 6379)

len(pairs_actual) + len(pairs_synth)

7111

Check if there are NodeSets counted multiple times

len(pairs_synth)

6379

once = set()
twice = set()
for rs_hgs, rs_bnc in pairs_synth:
    assert rs_hgs.nodes is rs_bnc.nodes, (rs_hgs.nodes.name, rs_bnc.nodes.name)
    if rs_hgs.nodes.digest in once:
        if rs_hgs.nodes.digest in twice:
            print(f'more than twice: {rs_hgs.nodes.name}')
        twice.add(rs_hgs.nodes.digest)
    else:
        once.add(rs_hgs.nodes.digest)
len(once), len(twice)

(6379, 0)

Make DataFrame

schema_hgs = dict(
    id=int,
    synthetic=bool,
    T=int,
    capacity=int,
    detextra=float,
    length=float,
    # incumbent=float,
    # incumbent_bound=float,
    exact_bound=float,
    # incumbent_id=int,
)

data = []
for rs_hgs, rs_bnc in pairs_actual + pairs_synth:
    data.append(dict(
        id=rs_hgs.id,
        synthetic=rs_hgs.nodes.name[0] == '!',
        T=rs_hgs.T,
        capacity=rs_hgs.capacity,
        detextra=rs_hgs.detextra,
        length=rs_hgs.length,
        exact_bound=rs_bnc.misc['bound'],
    ))
    # print(details['incumbent_id'], end=' ')
df_hgs = pl.DataFrame(data, schema=schema_hgs)

df_hgs.columns

['id', 'synthetic', 'T', 'capacity', 'detextra', 'length', 'exact_bound']

df_hgs['synthetic'].value_counts()

shape: (2, 2)

synthetic	count
bool	u32
true	6379
false	732

df_hgs['detextra'].describe()

shape: (9, 2)

statistic	value
str	f64
"count"	5129.0
"null_count"	1982.0
"mean"	0.003079
"std"	0.005495
"min"	-0.018452
"25%"	0.000774
"50%"	0.001803
"75%"	0.003515
"max"	0.118956

detextra_hgs = np.sort(np.nan_to_num(np.clip(df_hgs['detextra'], 0., None), nan=0.))

# Negative detextra in hgs is OK.
# It means some subtree roots chosen by PathFinder were closer to root
# than the path's head (since PathFinder relaxes the path-subtree constraint)
df_hgs.filter(pl.col('detextra') < -1e-3)

shape: (35, 7)

id	synthetic	T	capacity	detextra	length	exact_bound
i64	bool	i64	i64	f64	f64	f64
21006	false	60	10	-0.002975	51296.546654	49611.477292
21209	false	55	2	-0.003444	182175.265645	182804.81475
21208	false	55	4	-0.001266	99909.593562	99531.718807
21556	false	119	2	-0.017189	405793.043417	412587.411775
21557	false	119	3	-0.007867	291621.821448	292531.694887
…	…	…	…	…	…	…
16186	true	70	10	-0.009745	7.983239	7.73576
16857	true	50	4	-0.003896	8.106149	8.051141
17139	true	133	10	-0.002333	12.255932	11.63741
17696	true	68	8	-0.012723	8.916369	8.541159
17764	true	57	4	-0.001816	10.979295	10.874424

B&C solver

Make DataFrame

schema_exact = dict(
    id=int,
    synthetic=bool,
    T=int,
    capacity=int,
    detextra=float,
    length=float,
    bound=float,
)

data = []
for rs_hgs, rs_bnc in pairs_actual + pairs_synth:
    data.append(dict(
        id=rs_bnc.id,
        synthetic=rs_bnc.nodes.name[0] == '!',
        T=rs_bnc.T,
        capacity=rs_bnc.capacity,
        detextra=rs_bnc.detextra,
        length=rs_bnc.length,
        bound=rs_bnc.misc['bound'],
    ))
df_BnC = pl.DataFrame(data, schema=schema_exact)

df_BnC.height

7111

df_BnC['detextra'].describe()

shape: (9, 2)

statistic	value
str	f64
"count"	4935.0
"null_count"	2176.0
"mean"	0.002224
"std"	0.003292
"min"	-0.018452
"25%"	0.000616
"50%"	0.001563
"75%"	0.002992
"max"	0.120802

detextra_BnC = np.sort(np.nan_to_num(np.clip(df_BnC['detextra'], 0., None), nan=0.))

3.3 Impact of detours on the solver’s total length

with plt.style.context('pub_paperfarm'):
    data = 100*detextra_BnC

    fig, (ax1, ax2) = plt.subplots(2, 1, height_ratios=[1, 5], figsize=(3.5, 2),
                                   facecolor='w', sharex=True)
    # fig.subplots_adjust(hspace=0.05)
    ax1.spines.bottom.set_visible(False)
    ax1.xaxis.tick_top()
    ax1.tick_params(labeltop=False)
    num_data = len(data)
    below1 = data <= 1
    below2 = data <= 2
    bins = np.linspace(0, 3, num=76)
    # bins = None

    # top y_axis segment
    H1, X1, art1_ = ax1.hist(data[below1], bins = bins, color='tab:blue',
                             histtype='stepfilled', alpha=0.3)
    H2, X2, art2_ = ax1.hist(data[below2], bins = bins, edgecolor='tab:blue',
                             histtype='stepfilled', hatch=r'\\\\', facecolor='none')
    H, X, art_ = ax1.hist(data, bins = bins, color='tab:blue', histtype='step')

    # main y_axis segment
    H1, X1, art1_ = ax2.hist(data[below1], bins = bins, color='tab:blue',
                             histtype='stepfilled', alpha=0.3)
    H2, X2, art2_ = ax2.hist(data[below2], bins = bins, edgecolor='tab:blue',
                             histtype='stepfilled', hatch=r'\\\\', facecolor='none')
    H, X, art_ = ax2.hist(data, bins = bins, color='tab:blue', histtype='step')

    ax1.set_ylim(3100, 3200)  # outliers only
    ax2.set_ylim(0, 600)  # most of the data

    C1, C2 = H1.sum(), H2.sum()
    leg = ax2.legend(  # the empty Line2D is a hack to get a line instead of a patch
        [Line2D([], []), art2_[0], art1_[0]],
        ['Instance count',
         f'≤ 2%: {100*C2/num_data:.1f}%',
         f'≤ 1%: {100*C1/num_data:.1f}%'],
        title=f'B&C instances: {num_data}'
    )
    ax2.set_ylabel('Count')
    ax2.set_xlabel('Relative increase (%)')
    leg.legend_handles[0] = Line2D([], [])
    ax1.figure.savefig('fig_detextra_BnC.pdf', transparent=False)
    pass

with plt.style.context('pub_paperfarm'):
# with plt.style.context('pub_paperfarm'):
# if True:
    data = 100*detextra_hgs

    fig, (ax1, ax2) = plt.subplots(2, 1, height_ratios=[1, 5], figsize=(3.5, 2),
                                   facecolor='w', sharex=True)
    # fig.subplots_adjust(hspace=0.05)
    ax1.spines.bottom.set_visible(False)
    ax1.xaxis.tick_top()
    ax1.tick_params(labeltop=False)
    num_data = len(data)
    below1 = data <= 1
    below2 = data <= 2
    bins = np.linspace(0, 3, num=76)
    # bins = None

    # top y_axis segment
    H1, X1, art1_ = ax1.hist(data[below1], bins = bins, color='tab:blue',
                             histtype='stepfilled', alpha=0.3)
    H2, X2, art2_ = ax1.hist(data[below2], bins = bins, edgecolor='tab:blue',
                             histtype='stepfilled', hatch=r'\\\\', facecolor='none')
    H, X, art_ = ax1.hist(data, bins = bins, color='tab:blue', histtype='step')

    # main y_axis segment
    H1, X1, art1_ = ax2.hist(data[below1], bins = bins, color='tab:blue',
                             histtype='stepfilled', alpha=0.3)
    H2, X2, art2_ = ax2.hist(data[below2], bins = bins, edgecolor='tab:blue',
                             histtype='stepfilled', hatch=r'\\\\', facecolor='none')
    H, X, art_ = ax2.hist(data, bins = bins, color='tab:blue', histtype='step')

    ax1.set_ylim(2800, 2900)  # outliers only
    ax2.set_ylim(0, 600)  # most of the data

    C1, C2 = H1.sum(), H2.sum()
    leg = ax2.legend(  # the empty Line2D is a hack to get a line instead of a patch
        [Line2D([], []), art2_[0], art1_[0]],
        ['Instance count',
         f'≤ 2%: {100*C2/num_data:.1f}%',
         f'≤ 1%: {100*C1/num_data:.1f}%'],
        title=f'HGS instances: {num_data}'
    )
    ax2.set_ylabel('Count')
    ax2.set_xlabel('Relative increase (%)')
    leg.legend_handles[0] = Line2D([], [])
    ax1.figure.savefig('fig_detextra_HGS.pdf', transparent=False)
    pass

3.4 Comparison of B&C and HGS solutions

with plt.style.context('pub_paperfarm'):
    data_HGSvsBnC_gap = 1 - df_BnC['length']/df_hgs['length']
    data = 100*data_HGSvsBnC_gap.to_numpy()

    fig, ax = plt.subplots(facecolor='w', figsize=(3.5, 2))
    num_data = len(data)
    below1 = data <= 1
    below2 = data <= 2
    elsewhere = data > 2
    bins = np.linspace(-2, 8, num=201)

    # main y_axis segment
    H1, X1, art1_ = ax.hist(data[below1], bins = bins, color='tab:blue',
                             histtype='stepfilled', alpha=0.3)
    H2, X2, art2_ = ax.hist(data[below2], bins = bins, edgecolor='tab:blue',
                             histtype='stepfilled', hatch=r'\\\\', facecolor='none')
    H, X, art_ = ax.hist(data, bins = bins, color='tab:blue', histtype='step')

    C1, C2 = H1.sum(), H2.sum()
    leg = ax.legend(  # the empty Line2D is a hack to get a line instead of a patch
        [Line2D([], []), art2_[0], art1_[0]],
        ['HGS to B&C gap',
         f'≤ 2%: {100*C2/num_data:.1f}%',
         f'≤ 1%: {100*C1/num_data:.1f}%'],
        title=f'Instances: {num_data}'
    )
    ax.set_ylabel('Count')
    ax.set_xlabel('Relative difference (%)')
    leg.legend_handles[0] = Line2D([], [])
    ax.set_xlim(-2, 8)
    ax.figure.savefig('fig_HGS_to_BnC_gap.pdf', transparent=False)
    pass

3.5 Dealing with crossings from HGS-CVRP

Select HGS-CVRP entries

df_hgs.height

7111

method_hgs = db.Method.select(
    lambda m:
    m.funname == 'hgs_cvrp'
    and m.options['timeLimit'] == 120
)

method_hgs.count()

2

hgs_all = db.RouteSet.select(
    lambda rs:
    rs.method in method_hgs
)
hgs_all.count()

7044

hgs_invalid = hgs_all.filter(
    lambda rs:
    rs.valid is False
)
hgs_invalid.count()

0

# just the instances repaired after the first hgs run
hgs_repaired = hgs_all.filter(
    lambda rs:
    rs.valid is None  # provisory valid marker
    and 'repaired' in rs.misc
    and 'hgs_reruns' not in rs.misc
)
hgs_repaired.count()

332

hgs_reruns = hgs_all.filter(
    lambda rs:
    rs.valid is None  # provisory valid marker
    and 'hgs_reruns' in rs.misc
)
hgs_reruns.count()

98

Statistics

hgs_all = db.RouteSet.select(
    lambda rs:
    rs.creator == 'baselines.hgs'
    # rs.method.solver_name == 'HGS-CVRP'
    # and rs.method.options['timeLimit'] == 120
)
hgs_all.count()

7111

# just the instances repaired after the first hgs run
hgs_repaired = hgs_all.filter(
    lambda rs:
    rs.valid is None  # provisory valid marker
    and 'repaired' in rs.misc
    and 'hgs_reruns' not in rs.misc
)
hgs_repaired.count()

332

hgs_reruns = hgs_all.filter(
    lambda rs:
    rs.valid is None  # provisory valid marker
    and 'hgs_reruns' in rs.misc
)
hgs_reruns.count()

98

(hgs_repaired.count() + hgs_reruns.count())/hgs_all.count()

0.06046969483898186

hgs_repaired.count()/hgs_all.count()

0.04668822950358599

hgs_reruns.count()/hgs_all.count()

0.013781465335395865

Examining re-runs

subcounts = []
for reruns in range(1, 10):
    count = hgs_reruns.filter(
        lambda rs:
        rs.misc['hgs_reruns'] == reruns
    ).count()
    subcounts.append(count)
    print(reruns, count, count/hgs_all.count())

1 81 0.011390802981296582
2 12 0.0016875263675994938
3 2 0.0002812543945999156
4 1 0.0001406271972999578
5 2 0.0002812543945999156
6 0 0.0
7 0 0.0
8 0 0.0
9 0 0.0

total = sum(subcounts)
for reruns, count in zip(range(1, 10), subcounts):
    print(reruns, f'{100*count/total:.1f}%')

1 82.7%
2 12.2%
3 2.0%
4 1.0%
5 2.0%
6 0.0%
7 0.0%
8 0.0%
9 0.0%