MIP ile Çözüm — PuLP ve CBC

Geçen yazı bir SQL sorgusu ve on iki satırlık bir döngüyle çözücünün bulduğu değerin yüzde 99,4’ünü buldu. Şimdi kalan yüzde 0,6 için ne gerektiğine bakacağız: altmış yedi satır Python ve bedava bir çözücü.

Yalnız baştan söyleyeyim, bu yazının asıl konusu o yüzde 0,6 değil. Geçen yazının son cümlesi şuydu: yüzde 99,4’ü 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ı.

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

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]
    model += (
        pulp.lpSum(int(adaylar.loc[i, "adet"]) * x[i] for i in grup.index)
        >= p.min_koli * y[r]
    )

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 “Trabzon’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ı (5.006.978,84 TL) çıkarken her koşuda planın yirmi civarında hareketi farklıydı. Aday tablosunu sıralamak sorunu kapattı. Ders şu: determinizm tek bir yerde sağlanmaz, veriden çözüme kadar bütün zincirde sağlanır — ve bir yerde kaçırdığınızda size hata olarak değil, “geçen hafta başka söylüyordu” olarak geri döner.

Çözdürmek ve gelen cevabı okumak

sonuc = model.solve(pulp.PULP_CBC_CMD(msg=0, timeLimit=p.mip_zaman_limiti_sn))
if sonuc == pulp.LpStatusOptimal:
    durum = "optimal"
elif sonuc == pulp.LpStatusNotSolved:
    durum = "limit"
else:
    durum = "hata"

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

• optimal — bulunan çözümün daha iyisi yok, kanıtlanmış.
• limit — süre doldu. Elde bir çözüm var ve muhtemelen iyi, ama daha iyisi olmadığı kanıtlanmadı.
• hata — çözücü bir cevap üretemedi.

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.

MIP tam olarak nerede kazanıyor

Geçen yazıda iki planı yan yana koymuş ve şunu bulmuştuk: ikisi de aynı 252 option’ı, aynı 4.785 adedi taşıyor; fark yalnız rota sayısında — 196’ya karşı 140, aradaki 56 rota × 500 TL = 28.000 TL. Şimdi o farkın kodda nereden çıktığına bakalım.

Rotaların kaç option taşıdığına bakınca fark görünür hale geliyor:

Açgözlü planda 147 sevkiyat, tek bir option taşımak için açılıyor. Çözücüde bu sayı 56. Çözücü aynı işi daha az kamyonla yapı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 rotası:

M025 → M022 · 4 option · 89 adet

  MDL005-LCV   açgözlü planda:  M025 → M004
  MDL006-SYH   açgözlü planda:  M025 → M003
  MDL007-BEJ   açgözlü planda:  M025 → M022
  MDL071-BEJ   açgözlü planda:  M025 → M010

Çözücü bu dört bloğu tek kamyona koyuyor; açgözlü yöntem aynı dört bloğu dört ayrı sevkiyata dağıtmış. Dikkat: açgözlü yöntem yanlış bir şey yapmıyor — her bloğu tek tek en yüksek puanlı hedefine gönderiyor. Dört bloğun her biri kendi başına daha iyi eşleşmiş olabilir; ama dördü için dört kez kamyon parası ödeniyor ve toplamda geride kalıyor.

Bunun ikinci bir sonucu var, tabloda görünmüyor: iki plan aynı 252 bloğu kaldırıyor — aynı mağazadan, aynı option — ama bunların 111’inde alıcı farklı. Yani çözücü aynı planı ucuzlatmakla kalmıyor, kimin kime göndereceğini de yeniden kuruyor.

Neden yarım saniyeden az?

Üçüncü yazıda bu modelin 1.058 ikili değişkeni olduğunu ve kaba kuvvetle $2^{1058}$ olasılık demek olduğunu yazmıştım. Çözücünün bunu yarım saniyeye varmadan 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:

Boşluk sıfır. İlk tepkim “ölçüm bozuk” oldu, o yüzden gevşetmenin gerçekten gevşediğine ayrıca baktım: gevşetilmiş çözümde 664 $x$ değişkeninin 12’si ve 394 $y$ değişkeninin 6’sı gerçekten kesirli çıkıyor. Yani çözücü tam sayı olmayan bir tepe noktasında duruyor; buna rağmen o noktanın değeri tam sayı optimumla kuruşu kuruşuna aynı.

Anlamı şu: çözücü daha ilk adımda “en iyi çözüm 5.006.978,84’ten fazla olamaz” diye bir tavan çiziyor ve elinde tam o tavana oturan bir tam sayı çözüm oluyor. Tavanın aşağı inmesi hiç gerekmiyor. Kaç dal açıldığını ölçmedim, o yüzden “ağaç hiç açılmadı” demiyorum; bildiğim tek şey, çözücü işini bitirdiğinde kapatılacak boşluğun kalmadığı. Ama bunu “işte dal-sınırın kanıtı” diye satmak da yanlış olur. Bu sayı yöntemin gücünü değil, bu haftanın fotoğrafının uysallığını ölçüyor. Boşluğun sıfır çıkması bir garanti değil; başka bir hafta, başka bir eşik takımı aynı modelde o boşluğu pekâlâ açabilir. Dürüst cümle şu: bu örnekte çözücünün budayacak bir şeyi yoktu.

Her MIP böyle uslu değildir

Buradan “çözücüler her zaman hızlıdır” sonucunu çıkarmayın; bu problem kolay bir MIP. Kolaylığın sebepleri belli: LP boşluğu kapanmış durumda, kısıtlar birbirini az kilitliyor, ve aday kümesi zaten 664 satıra indirilmiş durumda.

Aynı model şu üç durumda çok daha zorlanır: aday kümesi elenmeden verilirse (144.000 satır), kapasiteler daralıp kısıtlar birbirini kilitlemeye başlarsa, ya da rota sabiti büyüyüp “hangi rotalar açılsın” sorusu asıl karar hâline gelirse. O zaman zaman limiti gerçekten devreye girer ve limit bayrağıyla gelen çözümle çalışmayı öğrenmek gerekir: çözüm iyidir, kanıtlanmamıştır ve raporda böyle söylenmelidir.

Bedelin dürüst muhasebesi

PuLP bedava, CBC bedava, çalıştırmak ~0,4 saniye. O zaman bu tartışma neden var?

Çünkü asıl maliyet lisansta değil. Modeli kuran kişi gittiğinde geriye kalan şey altmış yedi 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ığı, 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 arasında toplamı en büyük olandı.” 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 daha çok para kazandıran plana geçtiğinizde, boş çıkan müşterinin daha azını kurtarıyorsunuz 🙂