TSPLIB verisini MATLAB ile işlenebilecek hale nasıl getirebiliriz?

Sayın Mustafa Servet KIRAN hocamın aşağıdaki yorumu sonucu yazı içeriği değişmiştir.

EDGE_WEIGHT_TYPE: GEO
alanı oldukça önemli, GEO ise farkı denklem ile, EUC ise öklit denklemi ile uzaklık hesaplanmalıdır. GEO’yu EUC şeklinde hesaplarsanız “optimumdan daha optimum” (!) sonuç elde edebilirsiniz.
Cevahir her zamanki gibi detaylara önem vermeden ana mantık faydalı olmaktır prensibiyle yukarıdaki paylaşımı yaptığı için böyle bir hatanın oluşması kaçınılmazdı.

Bu arada pdist fonksiyonu matlabda çiftler arasındaki öklit mesafesini hesaplar, yukarıdaki hesaplama dolayısıyla yanlıştır.

Gezgin satıcı problemine çözüm ararken, literatürde bulunan bir çok kıyas problemi üzerinde test yapılmaktadır. TSPLIB tarafından sağlanan formatın kendi çözümlerimde kullanmak üzere aşağıdaki şekilde hazırlıyorum.

Öncelikle verilen koordinatların hangi tipte olduğuna dikkat etmemiz gerekiyor. Ben genelde öklid uzayında verilmiş verilerle çalıştığımdan, yukarıdaki uyarı yorumunu o yüzden aldım. Şimdi hem öklid hem de geo nasıl hesaplanır birlikte inceleyelim.

Aşağıdaki problem (ulysses22.tsp), EDGE_WEIGHT_TYPE: GEO olduğu için, verilen koordinatlar enlem ve boylam bilgisini içermektedir.

NAME: ulysses22.tsp
TYPE: TSP
COMMENT: Odyssey of Ulysses (Groetschel/Padberg)
DIMENSION: 22
EDGE_WEIGHT_TYPE: GEO
DISPLAY_DATA_TYPE: COORD_DISPLAY
NODE_COORD_SECTION
 1 38.24 20.42
 2 39.57 26.15
 3 40.56 25.32
 4 36.26 23.12
 5 33.48 10.54
 6 37.56 12.19
 7 38.42 13.11
 8 37.52 20.44
 9 41.23 9.10
 10 41.17 13.05
 11 36.08 -5.21
 12 38.47 15.13
 13 38.15 15.35
 14 37.51 15.17
 15 35.49 14.32
 16 39.36 19.56
 17 38.09 24.36
 18 36.09 23.00
 19 40.44 13.57
 20 40.33 14.15
 21 40.37 14.23
 22 37.57 22.56
EOF

TSPLIB kendi sayfasında bu problem için optimum değeri 7013 olarak vermiş, yani bizim hesaplamamız optimum tur değerini aldığı zaman bu değeri vermelidir.

function [ toplamuzaklik ] = objgeo()
latitude=[38.2400000000000;39.5700000000000;40.5600000000000;36.2600000000000;33.4800000000000;37.5600000000000;38.4200000000000;37.5200000000000;41.2300000000000;41.1700000000000;36.0800000000000;38.4700000000000;38.1500000000000;37.5100000000000;35.4900000000000;39.3600000000000;38.0900000000000;36.0900000000000;40.4400000000000;40.3300000000000;40.3700000000000;37.5700000000000];
longitude=[20.4200000000000;26.1500000000000;25.3200000000000;23.1200000000000;10.5400000000000;12.1900000000000;13.1100000000000;20.4400000000000;9.10000000000000;13.0500000000000;-5.21000000000000;15.1300000000000;15.3500000000000;15.1700000000000;14.3200000000000;19.5600000000000;24.3600000000000;23;13.5700000000000;14.1500000000000;14.2300000000000;22.5600000000000];
latitude = deg2rad(latitude);
longitude = deg2rad(longitude);
D=22;
birey=[1 14 13 12 7 6 15 5 11 9 10 19 20 21 16 3 2 17 22 4 18 8 1];
toplamuzaklik=0;
RRR = 6378.388;
for i=1:D  
    q1 = cos(longitude(birey(i)) - longitude(birey(i+1)));
    q2 = cos(latitude(birey(i)) - latitude(birey(i+1)));
    q3 = cos(latitude(birey(i)) + latitude(birey(i+1)));
    dij = round((RRR * acos( 0.5*((1.0+q1)*q2 - (1.0-q1)*q3) ) + 1.0));   
    toplamuzaklik=toplamuzaklik+dij;
end
end

kodunu çalıştırdığımız zaman 6971 değerini vermektedir.

Öncelikle yukarıdaki hesaplamayı TSPLIB’ın kendi sayfasından aldım.

Farklı olan derecenin radyana dönüşümü. O dönüşümden dolayı bir fark olabilir mi diye ilgili sayfadaki derece radyan dönüşümünü yaptım.

function [y] = degree2radian(x)
[~,D]=size(x');
PI = 3.141592;
for i=1:D   
    deg = round(x(i));
    min = x(i)- deg;
    rad(i) = PI * (deg + 5.0 * min/ 3.0) / 180.0;    
end
y=rad;
end

Bu dönüşüm ile yaptığım hesaplamada ise 7132 buldum.

Dikkat ederseniz RRR = 6378.388 şeklinde bir sabit, farklı kaynaklarda bu sayının farklı alındığını gördüm. O yüzden bir hesaplama farklılığı olabilir.

Şimdi gelelim öklid uzayında verilmiş bir problemin dönüştürülmesine:

eil51 problemi aşağıdaki gibidir:

NAME : eil51
COMMENT : 51-city problem (Christofides/Eilon)
TYPE : TSP
DIMENSION : 51
EDGE_WEIGHT_TYPE : EUC_2D
NODE_COORD_SECTION
1 37 52
2 49 49
3 52 64
4 20 26
5 40 30
6 21 47
7 17 63
8 31 62
9 52 33
10 51 21
11 42 41
12 31 32
13 5 25
14 12 42
15 36 16
16 52 41
17 27 23
18 17 33
19 13 13
20 57 58
21 62 42
22 42 57
23 16 57
24 8 52
25 7 38
26 27 68
27 30 48
28 43 67
29 58 48
30 58 27
31 37 69
32 38 46
33 46 10
34 61 33
35 62 63
36 63 69
37 32 22
38 45 35
39 59 15
40 5 6
41 10 17
42 21 10
43 5 64
44 30 15
45 39 10
46 32 39
47 25 32
48 25 55
49 48 28
50 56 37
51 30 40

Yukarıdaki koordinat bilgilerini Matlab içerisine aktarıyoruz.
a=[51×2 boyutunda koordinatların bulunduğu matris];
b=pdist(a); / Koordinatların birbirine göre uzaklıkları hesaplanıyor.
c=squareform(b); / Uzaklıkları 51×51’lik kare matris haline çeviriyoruz.

optimum 426 iken benim hesabımda 429.9833 bulunuyor. Sebebi ondalıklı sayılardan kaynaklı hassasiyettir.

Bu arada pdist fonksiyonu matlabda çiftler arasındaki öklit mesafesini hesaplar, yukarıdaki hesaplama dolayısıyla yanlıştır.

pdist komutu 2 sütun şeklinde verilen verileri de hesaplıyor.