MIP ile Çözüm — PuLP ve CBC

Geçen yazı bir SQL sorgusu ve kısa bir döngüyle çözücünün bulduğu değerin %90,9’unu buldu. Şimdi kalan %9,1 için ne gerektiğine bakacağız: kırk satıra yakın Python ve bedava bir çözücü.

Yalnız baştan söyleyeyim, bu yazının asıl konusu o %9,1 değil. Geçen yazının son cümlelerinden biri şuydu: %90,9’u biz biliyoruz, çünkü diğerini de kurup ölçtük. Modeli kurmanın esas getirisi bu: kaybınızın büyüklüğünü öğreniyorsunuz. Kâr olarak geri dönen kısım cabası; gerçi bu hafta o caba da küçük değil, 3.253.435 TL.

PuLP: matematiği Python’da yazmak

Çözücüler modeli standart bir biçimde ister; kimse o biçimi elle yazmaz. Arada bir modelleme dili durur. Python tarafında en yaygını PuLP:

pip install pulp

Bu kadar. Çözücü de (CBC, COIN-OR’un açık kaynak dal-sınır motoru) kutunun içinde geliyor; ayrıca lisans, sunucu, kurulum yok. Modelin tamamı aşağıdaki gibi kuruluyor ve her parçasının üstüne üçüncü yazıdaki formülü koyuyorum ki çeviriyi görebilesiniz.

Değişkenler. $x_{ijo} \in \{0, 1\}$ her aday hareket için, $y_{ij} \in \{0, 1\}$ her rota için:

model = pulp.LpProblem("blok_transfer", pulp.LpMaximize)
x = {i: pulp.LpVariable(f"x_{i}", cat="Binary") for i in adaylar.index}
rotalar = sorted(set(zip(adaylar.verici, adaylar.alici)))
y = {r: pulp.LpVariable(f"y_{r[0]}_{r[1]}", cat="Binary") for r in rotalar}

Amaç fonksiyonu. $\max \sum w_{ijo}x_{ijo} - R\sum y_{ij}$:

model += (
    pulp.lpSum(adaylar.loc[i, "w"] * x[i] for i in adaylar.index)
    - p.rota_sabiti_tl * pulp.lpSum(y.values())
)

Kısıt 1: blok tek hedefe. $\sum_j x_{ijo} \le 1$:

for (verici, option), grup in adaylar.groupby(["verici", "option_id"]):
    model += pulp.lpSum(x[i] for i in grup.index) <= 1          # blok tek hedefe

Kısıt 2: alıcı kapasitesi. $\sum s_{io}x_{ijo} \le Q_j$:

for alici, grup in adaylar.groupby("alici"):
    model += (
        pulp.lpSum(int(adaylar.loc[i, "adet"]) * x[i] for i in grup.index)
        <= kapasite.get(alici, 0)
    )

Kısıt 3 ve 4: rota açılışı ve minimum koli. $x_{ijo} \le y_{ij}$ ve $\sum_o s_{io}x_{ijo} \ge K y_{ij}$:

for r, grup in adaylar.groupby(["verici", "alici"]):
    for i in grup.index:
        model += x[i] <= y[r]                                    # rota açılmadan taşıma yok
    model += (
        pulp.lpSum(int(adaylar.loc[i, "adet"]) * x[i] for i in grup.index)
        >= p.min_koli * y[r]                                     # açık rota koliyi doldurur
    )

Hepsi bu. Formüllerle kod arasında birebir karşılık var; PuLP’ın yaptığı iş += ile eklediğiniz satırları toplayıp çözücünün anlayacağı biçime çevirmek.

Bir satırlık ders: sorted

Yukarıdaki değişken bloğunda küçük bir ayrıntı var:

rotalar = sorted(set(zip(adaylar.verici, adaylar.alici)))

O sorted başta yoktu. Kod incelemesinde eklendi ve sebebi şu: Python’da bir kümenin (set) dolaşma sırası çalışmadan çalışmaya değişebilir. Sıra değişince model dosyasındaki değişken sırası değişir; değişken sırası değişince, eşit değerli birden çok en iyi çözüm varsa çözücü bir seferinde birini, bir seferinde diğerini döndürür.

Matematiksel olarak ikisi de doğrudur; aynı parayı kazandırırlar. Ama pazartesi sabahı raporu alan Ali için doğru olmaları yetmez: geçen hafta “Malatya’dan al” diyen sistemin bu hafta aynı veriyle “Kayseri’den al” demesi, sistemin güvenilirliğini bitirir. Geçen yazıda açgözlü döngünün sıralamasına verici, alici, option_id kuyruğunu eklerken aynı derdi çözüyorduk. Determinizm bir matematik meselesi değil, güven meselesidir ve iki yöntemde de ayrı ayrı ele almak gerekiyor.

Bu yazıyı yazdıktan sonra aynı dersin ikinci örneği çıktı, o yüzden buraya ekliyorum: sorted tek başına yetmemiş. Aday tablosunu üreten SQL sorgularında order by yoktu ve DuckDB satırları her koşuda aynı sırada döndürmek zorunda değil. Sıra değişince değişkenlerin adı değişiyor, adlar değişince çözücü yine eşit değerli başka bir optimuma düşüyordu: art arda üç koşuda amaç değeri kuruşu kuruşuna aynı çıkarken her koşuda planın yirmi civarında hareketi farklıydı. Aday tablosunu sıralamak sorunu kapattı; geçen yazıdaki sorgunun son satırındaki order by oradan geliyor. Ders şu: determinizm tek bir yerde sağlanmaz, veriden çözüme kadar bütün zincirde sağlanır. Bir yerde kaçırdığınızda da size hata olarak değil, “geçen hafta başka söylüyordu” olarak geri döner.

Çözdürmek ve gelen cevabı okumak

model.solve(pulp.PULP_CBC_CMD(
    msg=0,
    threads=0,
    gapRel=p.mip_bosluk_orani,
    timeLimit=p.mip_zaman_limiti_sn,
    maxNodes=p.mip_dugum_limiti,
    logPath=log_yolu,
))

msg ve logPath çözücünün ne yazdığıyla ilgili; kalan dört ayarın her birinin bir hikâyesi var.

gapRel çözücüye ne zaman duracağını söylüyor: elindeki planla üst sınır arasındaki fark %0,5’in altına inince. Yani “en iyiyi bul” demiyoruz, “en iyiye yeterince yakın olduğunu kanıtladığın an dur” diyoruz.

maxNodes ile timeLimit emniyet sınırları ve aralarındaki fark yine determinizm dersi. Süre sınırı makinenin hızına bağlıdır: hızlı bilgisayar aynı sürede daha çok dal açar, başka bir plan bulur. Düğüm sınırı ise makineden bağımsızdır; o yüzden süreden önce bağlayacak şekilde ayarlı. Emniyete takılan bir model de, hangi bilgisayarda koşarsa koşsun, aynı planı verir. threads=0 da aynı sebeple orada: paralel arama aynı girdide kimi koşuda takılıyor, kimi koşuda çöküyordu; seri arama her seferinde aynı yolu yürüyor.

Bir şeyi de bilerek yapmadık: açgözlü planı çözücüye başlangıç noktası olarak vermedik. Denedik; CBC’nin bu sürümünde o yol bozuktu, düzeltince de sınıra takılan modelleri kötüleştirdi.

Çözücü durunca ne bulduğunu da ondan okumamız gerekiyor:

# PuLP'nin `solve()` dönüşü (LpStatus) süre ya da düğüm sınırında olurlu
# çözümle duran CBC için de LpStatusOptimal'dir: "optimal" kanıt demek
# değildir. Kanıtı çözüm durumu (`sol_status`) taşır.
DURUMLAR = {
    pulp.LpSolutionOptimal: "optimal",          # boşluk toleransı içinde (bkz. cozumle)
    pulp.LpSolutionIntegerFeasible: "limit",    # emniyet sınırında durdu, olurlu en iyi çözüm
}

durum = DURUMLAR.get(model.sol_status, "hata")

Buradaki üç yol, çözücüyü üretime koyarken en çok ihmal edilen yerdir. Çözücü her zaman “en iyisini buldum” demez:

• optimal: çözücü toleransın içinde durdu. Plan, en iyi plandan en fazla %0,5 uzak; ne kadar uzak olduğunu da logdaki son üst sınırdan okuyoruz. Kanıtlı en iyi demek değil.
• limit: bir emniyet sınırı doldu. Elde bir çözüm var ve muhtemelen iyi, ama tolerans kadar yakın olduğu bile kanıtlanmadı.
• hata: çözücü bir cevap üretemedi.

Bu tuzağa düşmek çok kolay: solve()’un dönüşüne bakan kod, süre sınırında yarıda kesilmiş planlara da sessizce “en iyi” der. Doğru bilgi sol_status’ta.

Açgözlü yöntemde böyle bir ayrım yok; o her zaman bir liste veriyor. Model kurmanın bedellerinden biri bu belirsizliği kodda karşılamak zorunda olmanız. Karşılamazsanız, günün birinde süre dolar ve elinizde sessizce yarım kalmış bir plan kalır; üstelik raporunuz ona “optimal” der.

MIP tam olarak nerede kazanıyor

Geçen yazıda farkın neredeyse tamamının mal seçiminden geldiğini gördük: 3.253.435 TL’nin yalnız 11.000 TL’si rotadan. Yine de rotalar işin görünen yüzü. Çözücü daha çok bloğu daha az sevkiyatla taşıyor ve rotaların kaç option taşıdığına bakınca bu görünür hâle geliyor:

Tabloda yalnız ilk dört satır var; kuyruk iki planda da uzun, en dolu rota ikisinde de 17 option taşıyor. Açgözlü planda 268 sevkiyat tek bir option taşımak için açılıyor. Çözücüde bu sayı 194. Çözücü daha çok malı daha az kamyonla taşıyor, çünkü amaç fonksiyonundaki $-R\sum y_{ij}$ terimi ona her yeni rotanın 500 TL’ye mal olduğunu söylüyor. Açgözlü döngüde böyle bir terim yok, olamaz da: döngü sıradaki hareketi görmediği için “bunu şuraya koysam kamyon zaten kalkacaktı” diye düşünemez.

Somut bir örnek, çözücünün en dolu rotalarından biri:

Kayseri 1 → İzmir 1 · 17 option · 5.722 adet

  MDL0540-BEJ   569 adet   açgözlü planda: Kayseri 1 → İstanbul 2
  MDL0541-KRM   236 adet   açgözlü planda: Kayseri 1 → Ankara 2
  MDL0543-HAK   478 adet   açgözlü planda: Kayseri 1 → İstanbul 2
  MDL0544-IND   372 adet   açgözlü planda: Kayseri 1 → İstanbul 2
  MDL0544-YSL   280 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0551-EKR   250 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0551-HAK   233 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0554-BRD   388 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0555-BRD   363 adet   açgözlü planda: Kayseri 1 → İstanbul 3
  MDL0556-KRM   202 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0557-SYH   607 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0576-BRD   127 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0576-KRM   189 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0582-IND   147 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0597-LCV   165 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0601-LCV   597 adet   açgözlü planda: Kayseri 1 → İzmir 1
  MDL0603-YSL   519 adet   açgözlü planda: Kayseri 1 → İstanbul 3

Çözücü bu on yedi bloğu tek kamyona koyuyor; açgözlü yöntem aynı blokları dört ayrı alıcıya dağıtmış: 11’i yine İzmir 1’e, kalanı İstanbul 2’ye, İstanbul 3’e ve Ankara 2’ye. Dikkat: açgözlü yöntem yanlış bir şey yapmıyor; her bloğu sırası geldiğinde kalan en iyi hedefe gönderiyor. Çözücünün o altı bloğu neden İzmir 1’e çektiğini ise tek bir rotaya bakarak okuyamazsınız. Rota masrafı bir sebep olabilir; İstanbul’daki iki mağazanın boşluğunu başka bloklara bırakmak başka bir sebep. Çözücü hiçbir hareketi tek başına seçmiyor, bu da yazının sonunda başımıza iş açacak.

Tabloda görünmeyen bir sonuç daha var: iki planın ortak taşıdığı 571 bloğun 193’ünde alıcı farklı; çözücü açgözlünün hiç dokunmadığı 90 bloğu da yerleştiriyor. Yani çözücü açgözlünün planını ucuza yeniden dizmiyor; hangi malın taşınacağını da, kimin kime göndereceğini de baştan kuruyor.

Çözücü nerede duruyor?

Üçüncü yazıda bu modelin 7.730 ikili değişkeni olduğunu ve kaba kuvvetle $2^{7.730}$ olasılık demek olduğunu yazmıştım. Çözücünün bunu bir dakikanın altında bitirmesinin sebebini de anlatmıştım: dal-sınır. Şimdi ölçelim.

Modeli iki kez çözelim. Önce değişkenlerin 0 ile 1 arasında kesirli olmasına izin vererek (LP gevşetmesi), sonra gerçek hâliyle:

Gevşetmeyi çözmek tam sayı modele göre çok kısa sürüyor. Çıkan sayı bir tavan: hiçbir tam sayı plan 36.255.104,44 TL’yi geçemez. Gevşetilmiş çözümde 5.369 $x$ değişkeninin 105’i, 2.361 $y$ değişkeninin 56’sı kesirli çıkıyor; çözücü “bu bloğun yarısını gönder, bu rotayı yarım aç” diyor. Tavanla elimizdeki plan arasında %0,89’luk bir boşluk var; sürenin geri kalanı o boşluğu kapatmaya gidiyor.

Çözücü dallanarak tavanı aşağı çekiyor. Durduğunda son üst sınırı 35.997.742,47 TL: LP’nin bıraktığı boşluğun 257.361,97 TL’sini kapatmış. Elindeki plan bu sınırın %0,18 altında; tolerans %0,5 olduğu için orada duruyor. Bu yüzden bu yazıda hiçbir yerde “en iyi plan” demiyorum. Dürüst cümle şu: bu plan, olası en iyi plandan en fazla %0,18 uzak. Belki tam da en iyisi; bunu kanıtlamak için çözücüyü daha uzun koşturmak gerekirdi ve en fazla %0,18’lik bir iyileşme için o süreyi harcamadık.

Her MIP uslu değildir

Buradan “çözücüler her zaman hızlıdır” sonucunu çıkarmayın. Bu model referans ayarında toleransın içinde durdu ama kadranı çevirince öyle olmuyor. Dizinin sonunda karşılaştıracağımız altı eşik takımından dördünde çözücü emniyet için koyduğumuz düğüm sınırına takıldı ve limit bayrağıyla döndü. O planlar optimumdan en fazla %1,33 ile %5,24 arasında uzak; hepsi yine de açgözlüden fazla kazandırıyor, ama gerçek en iyinin o aralığın neresinde olduğunu bilmiyoruz.

Bir MIP’i neyin zorlaştırdığı aşağı yukarı bilinir; bunlardan ikisi bu veride var. Aday kümesi büyüyor: alıcı tavanı açıldıkça kırık olmayan raflar da alıcı olabiliyor. Kapasiteler bağlıyor: 80 mağazanın 12’sinin hiç boşluğu yok; referans planın dışında kalan puanı pozitif bloklardan 194’ü alıcısının kalan boşluğuna sığmıyor. Kısıtlar birbirini kilitledikçe tavan yavaş iniyor.

Böyle bir durumda yapılacak şey çözümü atmak değil, limit bayrağıyla çalışmayı öğrenmek: çözüm iyidir, kanıtlanmamıştır ve raporda böyle söylenmelidir. Bizim raporumuz her MIP satırının yanına durumunu ve boşluğunu yazıyor. Sizin verinizde zor olan büyük ihtimalle burası olacak: çözücünün söylediği “limit”i, planı kullanacak kişiye “yanlış” diye değil “kanıtsız” diye anlatmak.

Bedelin dürüst muhasebesi

PuLP bedava, CBC bedava, çalıştırmak bir dakikadan kısa. Açgözlünün 39,7 ms’sinin bin katından fazla, ama bir pazartesi sabahı için hiçbir şey. O zaman bu tartışma neden var?

Çünkü asıl maliyet lisansta da sürede de değil. Modeli kuran kişi gittiğinde geriye kalan şey kırk satır Python değil, o satırların neden öyle olduğunu bilen kimsenin kalmaması. Kapasite tanımı değiştiğinde kısıtı kimin güncelleyeceği, çözücü sürümü yükselince sonuçların neden bir tık kaydığı, limit dönen bir haftada planın kullanılıp kullanılmayacağı ve hepsinden zoru: mağaza müdürü “benim malımı niye aldınız?” diye sorduğunda verilecek cevap.

Açgözlü yöntemde o cevap bir satır: “puanı en yüksek eşleşme buydu.” Çözücüde cevap dürüstçe şudur: “seçilen plan, kısıtları sağlayan bütün planların en iyisinden en fazla %0,18 geride kalan plandı.” Doğru bir cevap ama ikna edici değil, çünkü tek bir hareketin gerekçesi diye bir şey yok; hareketler birbirine bağlı bir bütün olarak seçiliyor. Bu, optimizasyon projelerinin teknik değil, güven sebebiyle raflarda kalmasının en yaygın sebebi.

Sonraki yazıda iki yöntemi bütün ölçütlerle yan yana koyacağız: para, süre, kurulum, açıklanabilirlik. Orada beklemediğim bir sonuç çıktı: kadranı çevirip parayı üçte iki büyüten plana geçtiğinizde, boş çıkan müşterinin kurtarılan payı neredeyse yerinde sayıyor 🙂