sample.vehicle_routing のソースコード

# -*- coding: utf-8 -*-
from itertools import product
from math import ceil, sqrt
#from typing import override  # Python 3.12 or more
from random import sample, seed; seed(0)

from pysimple import Element, Parameter, IntegerVariable, BinaryVariable, Problem, NuoptStatus, Sum
from pysimple.typing import Table
from .column_generator import ColumnGenerator

class VehicleRouting(ColumnGenerator):
    """Vehicle Routing Problem"""
    def __init__(self, capacity: int, volume: dict[str, int], dis: dict[tuple[str, str], float]) -> None:
        self.capacity = capacity
        self.volume = volume  # volume[i]
        self.dis = dis        # dis[ijX]

    def create_init_pattern(self) -> None:
        """倉庫 X と店舗 i を往復するルート"""
        pattern, rd = {}, {}  # rd: ルート r の巡回距離
        for r, i in enumerate(self.volume, 1):
            pattern[r,i] = 1
            rd[r] = self.dis['X',i]*2
        super().create_init_pattern(pattern=pattern, c=rd)

    def create_new_pattern(self, lmbval: Table, *, silent: bool=True) -> tuple[Table, Table]:  # CVRP
        """Sum(pattern[p,i]*lmb[i], i) < rd[p] となる p を見つける"""
        iX = Element(value=['X', *self.volume])   # 倉庫(i=X)および店舗
        jX = Element(set=iX.set)
        ijX = iX!=jX
        i = iX!='X'  # 倉庫以外の店舗
        j = Element(set=i.set)
        ij = i!=j
        v = Parameter(index=i, value=self.volume)  # 店舗 i への配送量
        d = Parameter(index=ijX, value=self.dis)
        lmb = Parameter(index=j, value=lmbval)
        N = len(i.set)

        z = BinaryVariable(index=ijX)  # iX から jX に移動するか
        u = IntegerVariable(index=i, lb=1, ub=N)
        prb = Problem(type=max)
        prb += Sum(lmb[ij(1)]*z[ij]) - Sum(d[ijX]*z[ijX]), 'obj'
        prb += Sum(z['X',i], i) == 1,  '倉庫から出発'
        prb += Sum(z[i,'X'], i) == 1,  '倉庫へ帰還'
        prb += Sum(z[ij], ij(1)) <= 1, '店舗 i から出発は 1 台まで'
        prb += Sum(z[ij], ij(0)) <= 1, '店舗 j へ到着は 1 台まで'
        prb += Sum(z[ijX], ijX(0))[i] == Sum(z[ijX], ijX(1))[i], '店舗 i からの到着==出発'
        prb += u[ij(0)] - u[ij(1)] + N*z[ij] <= N-1, 'MTZ'
        iXj = iX!=j
        prb += Sum(v[iXj(1)]*z[iXj]) <= self.capacity, '最大積載量'
        prb.solve(silent=silent)
        assert prb.status == NuoptStatus.OPTIMAL
        return z[z[ijX].val>0].val, prb.objective.val

    #@override
    def update_pattern(self, zval: Table) -> None:
        """update self.pattern, self.c"""
        zdct: dict[str, str] = dict(zval.keys())
        route = [cur:='X'] + [cur:=zdct[cur] for i in range(len(zdct))]
        pn = self.patternnum
        print(f'new pattern{pn}', '->'.join(route))
        self.pattern.update({(pn, i): 1 for i in route[1:-1]})
        self.c[pn] = sum(self.dis[i,j] for i, j in zip(route[:-1], route[1:]))

    def visualize(self, zval: Table) -> None:
        """結果表示"""
        for r, in zval:
            route = [i for r_, i in self.pattern if r_==r]
            print(f'route{r}:', '->'.join(['X', *route, 'X']))

def create_init_data(*, N: int) -> tuple[int, dict[str, int], dict[tuple[str, str], float]]:
    """N: 店舗数"""
    mapsize = 2*ceil(sqrt(N))  # 密度 25% 程度
    I = [f'p{t}' for t in range(1, N+1)]  # 店舗の種類
    volume = dict(zip(I, sample(range(1, N+1), N)))  # 配送量 >0
    capacity = sum(volume.values())//2  # sum(volme) 以上だと一回で巡回できてしまう
    # 倉庫 X はマップの中央，他店舗はランダム
    point: list[tuple[float, ...]] = sample(list(product(range(mapsize+1), repeat=2)), N)
    I.insert(0, 'X')
    point.insert(0, (mapsize//2+0.5, mapsize//2+0.5))
    dis = {(i1,i2): sqrt((x1-x2)**2 + (y1-y2)**2)
           for (i1,(x1,y1)), (i2,(x2,y2)) in product(zip(I, point), repeat=2) if i1 != i2}
    return capacity, volume, dis


if __name__ == '__main__':
    CAPACITY, VOLUME, DIS = create_init_data(N=15)
    vr = VehicleRouting(CAPACITY, VOLUME, DIS)
    ZVAL = vr.solve()
    vr.visualize(ZVAL)