The Solution of Traveling Salesman Problem with Time Windows by
MILP in Matlab Code
JAROMÍR ZAHRÁDKA
Department of East Bohemia Region,
Žampach 10, 53401 Žamberk,
CZECH REPUBLIC
Abstract: - A delivery driver uses his truck to distribute ordered goods from the warehouse (depot) to
customers. Each customer determines the delivery point (by GPS coordinates) and time window for the
distribution of goods. This problem can be called the traveling salesman problem with time windows (TSPTW).
The objective of the solution is to select the sequence of delivery points so that the travel distance and the total
travel time are minimal. The necessary condition is that the deliverer leaves the warehouse, visits all delivery
points (where he arrives at the required time windows), and returns to the warehouse. In this article, one exact
robust solution of the TSPTW is presented by using mixed integer linear programming implemented in Matlab
code. The created algorithm can be used for any number of customers.
Key-Words: - Matlab, mixed-integer linear programming, MILP, optimization, traveling salesman problem
with time windows, TSPTW.
1 Introduction
This article aims to create a functional Matlab
algorithm that can find the optimal solution of the
traveling salesman problem with time windows
(TSPTW) for any given number of customers.
The simple traveling salesman problem (TSP)
lies in finding the shortest path or shortest time that a
salesman should take to visit all customers and
return to the depot. Related works on the topic of
TSP are [1], [2], [3], [4] and [5].
The TSPTW is an extension of the TSP which
requires, in addition, that the traveling salesman
arrive at the customers during pregiven time
windows.
Our work is focused on the use of branch and
bound, and branch and cut algorithms that have been
used, for example, in studies [5], [6] and [7]. The
TSPTW was solved in [8] and [9]. Authors used the
compressed-annealing heuristic and the QUBO
model.
The related problem of finding the real optimal
road route from Start-point to End-Point is described
by pseudocode for Dijkstra's algorithm in [10]. The
minimum distances of the routes are obtained using
google navigation which is also used in our work.
The TSP solution using Matlab code is
demonstrated in [11]. The Matlab code for the
solution of the traveling salesman problem was used
in [5] and [7].
2 Mixed-integer Linear Programming
The problem of mixed-integer linear programming
(MILP) is generally expressed by:
(1)
The vector is the column vector of all flow
variables. The symbol denotes the column vector
with coefficients represented in the objective
function which is . The number of components
of is equal to the number of all flow variables.
Writing means the list of variables of the
vector that takes only the integer values. The other
flow variables are treated as real variables.
The linear inequality constraint matrix is
specified as a matrix of real numbers, which are the
linear coefficients in the system of inequality
constraints given by . The size of the
matrix is , where is the number of
inequality constraints and is the number of flow
variables. Linear inequality constraints vector is
the vector of given real numbers (right sides of the
inequalities). The length of the vector is .
The linear equality constraint matrix is
specified as the matrix of real numbers, which are
the linear coefficients in the system of equations
. The size of the matrix is ,
Received: June 6, 2023. Revised: January 21, 2024. Accepted: February 8, 2024. Published: March 21, 2024.
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
186
Volume 21, 2024
where is the number of linear equation
constraints and is the number of flow variables.
Linear equality constraint vector is the vector of
given real numbers (right sides of the equations).
The length of the vector is .
The lower and upper bounds of flow variables
(items of ) are specified as real numbers - elements
of vectors , and . The vector satisfies the
inequality . The core of the solver for
the MILP solution in Matlab code is the command:
X=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub),
(2)
where previously used variables are expressed by
the corresponding variable identifiers of Matlab
langu-age. The transformation of all variable labels
into Matlab identifiers is contained in Table 1.
A more detailed explanation of the intlinprog
command is available in [11].
3 Mathematical Formulation of the
Traveling Salesman Problem with
Time Windows
The traveling salesman problem with time windows
(TSPTW) can be defined as follows. Let
be a connected graph consisting of a set of
nodes indexed by , and a set of non-
negatively weighted arcs between pairs of
corresponding nodes of the graph . The index
is used for the seller, and indexes for
customers. For easier reference, let be
the set of customers, and . Each of the
customers can be reached only within a
specified time interval (window) . Vectors
of lower and upper boundaries of the time window
are denoted by and . For each customer
let be the service time associated with the
handover and unloading of goods. The vector of
customer service times means . The
constant means the moment when the vehicle can
first leave the depot.
Let be the length of the one-way path (in
meters) and the one-way time-distance (in
seconds), i.e. the travel duration from node to node
for all . Therefore is a
distance matrix and is a time-distance
matrix. Both matrices and are non-negative,
asymmetric, and related (but essentially independent
non-negative), with zeros on the main diagonal, i.e.
and for each . It is necessary to
satisfy the triangular inequalities among all nodes of
the graph in both matrices.
Since we have to respect time windows, it is
necessary that the time-distance matrix be used
for optimization. The elements of the time-distances
matrix are not strictly proportional, in general, to
the corresponding elements of the distance matrix
. The driving time between two nodes is
significantly influenced by the quality of the roads
used.
To start the optimization process, it is necessary
to select the time moment when the vehicle leaves
the depot. Regarding the lower valid time windows
of all customers, the appropriate time is:
(3)
But when the first visited customer is
selected during the optimization process, for which
it is valid that , then a waiting time
will be required. To eliminate the
unnecessary waiting time at the first customer, it is
appropriate to shift the departure time of the vehicle
from the depot to the real value .
The core of the TSPTW solution is to find the
cycle in the graph which contains all nodes of
the graph so that the travel time is the shortest, and
where all time windows are respected. For this
purpose, integer variables for are
introduced, which can only take the values 0 or 1.
The value means that the arc from the node
to is included in the cycle. The value
means that the corresponding arc is not included.
For a systemic reason, variables are included,
but all these variables are fixed by the value zero
(), for each . Variables are elements
of a matrix . The number of all
variables is .
Other flow variables that need to be used in the
TSPTW solution are the vehicle departure times
from all customers. Therefore real (non-integer)
flow variables are introduced for each . The
variable indicates the moment when the driver
leaves the customer number . The variables are
arranged as elements of a vector .
So, the number of all flow variables is .
The solution of TSPTW is realized by the
optimal solution of the mixed-integer linear
programming problem:
subject to (4)
are integers, even true that (5)
(6)
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
187
Volume 21, 2024
(7)
(8)
(9)
(10)
(11)
(12)
In the above expressed model (4), with flow
variables and , the linear function:
(13)
is minimized. The first part,
guaran-
tees finding the cycle, which takes the minimum
travel time. For this, in the optimization process, it
is necessary that the second part of the minimized
function (13), i.e.
will be smaller by at
least 10 times. This condition is certainly met if we
choose the coefficients:
(14)
It is assumed that the values of the variables
do not exceed the order of 86400 sec (1 day).
This objective function guarantees that the
members containing the variables will not affect
the optimal values of the variables . At the same
time, the arrival times will be minimized and
downtimes of vehicles will not be necessary before
reaching customers.
Using statement (5) it is given that
are all integer variables, even binary ones.
If , the inequality (6) expresses the
relationship ,
for each . The right side of this
inequality can only be nonnegative because the
upper bound of time window plus service time
is greater than or equal to the departure time
(from the node ), and the departure time (from
the node ) is greater than or equal to the sum of
lower bound of the time window , and the service
time .
In the case the inequality (6) defines a
relationship . This
expresses the condition that the departure time
(from node ) is greater than or equal to the sum of
the departure time (from node ), the traveling
time from node to the one, and the service
time .
The Constraint (7) defines relations between
the flow variables and , for each . In the
case the inequality (7) expresses a
relationship , i.e. the departure time
from node is greater than or equal to the sum of
the lower bound of the time window and service
time . In the case the inequality expresses
the relationship . The departure
time from node is greater than or equal to the sum
of the departure time (from the depot), traveling
time (from the depot to the node ), and the
service time .
Statements (8), (9) and (10) declare
equation constraints. The inequalities in (11) declare
the lower and upper bounds (0 and 1) for variables
, , , . The inequalities in (12)
enforce permitted limits for departure times ,.
4 Transfer of TSPTW Model to the
Matlab Environment
The procedure for solving TSPTW in Matlab code is
contained in the M-function TSPTW_SOLVER.m,
which is listed as an Appendix at the end of this
article.
For solving the TSPTW in the Matlab
environment, it is necessary to reliably modify the
used variables following Matlab rules. The main
problem with the transformation is that the index 0
cannot be used in Matlab, therefore, some variables
need to be reindexed. An overview of the main
variables transformed into the Matlab is contained
in Table 1., where identifier p substitutes (n+1)^2.
We assume the vehicle goes around to
customers . Each customer requests the
arrival of a vehicle in a certain time window. The
lower and upper boundaries of the time windows
and service times for all customers are stored in M-
vectors lw, uw and m. For departure time from
the depot, the identifier t0 is used in Matlab. The
distance matrix and the travel time matrix have
to be transferred to the Matlab environment as
matrices D and C thus D(i,j) ≡ and
C(i,j) ≡ for indexes .
For the solution of TSPTW via the
intlinprog
command, all flow variables have to be arranged in
a column vector X. First flow variables of
vector X are elements of the matrix . Each
variable , is represented in Matlab code by
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
188
Volume 21, 2024
flow variable X(i*(n+1)+j+1,1). Each variable ,
is represented in Matlab code by flow variable
X((n+1)^2+i,1).
In the intlinprog command of Matlab the
objective function (13) is expressed like f'*X, where
f is a column vector with components.
The first components of the vector f are
elements of travel time matrix C such that f((i-
1)*(n+1)+j,1)=C(i,j), i,j , and the
last components have the same value according
to relation (14). The value of (kf in Matlab code)
is calculated by commands on lines No. 2 to 8 of the
Appendix. The vector f of all components of the
objective function are created in Matlab code on line
No. 9.
The value according to (3) is introduced in
Matlab code on row No. 10 of the Appendix.
The constraints (5), (6) are transformed to the
Matlab code by inequalities
(C(i+1,j+1)+uw(i)+m(i)-lw(i))*X(i*(n+1)+j+1,1)+…
X(p+i,1) - X(p+j,1)<= uw(i)-lw(j)+m(i)-m(j), for all
indexes i, j, i ≠ j, and
(C(1,j+1) +t0 -lw(j))*X(j+1,1) -X(p+j,1)<= -lw(j)-
m(j), for j .
Concerning the last two inequalities, we can fill
in the elements of matrix A and vector b, which are
essential input parameters of the intlinprog
function. The corresponding Matlab code can be
found on rows No. 11 to 24 of the Appendix.
Concerning equalities (7), (8), (9), we can fill the
elements into the matrix Aeq and vector beq. The
relevant Matlab code can be found on rows No. 25
to 43 of the Appendix.
Two input variables lb, ub of the intlinprog
command (2) are the lower and upper bounds of the
flow variables. Concerning the relations (9), (10),
(11) the relevant components are filled in using the
commands on lines No. 44 to 63 of the Appendix.
The vector of integer variables is, in the
Matlab, coded by the vector intcon, and specifies all
indices of flow variables, which are taken as
integers. It uses the command intcon=1:p on row
No. 64, p=(n+1)^2. This means that all flow
variables are taken as integers.
By installing the input variables f, intcon, A, b,
Aeq, beq, lb, and ub, and running the command
intlinprog (on row No. 65), we get the optimal
TSPWT solution if it exists. The optimal solution is
given by the optimal values of components of flow
variables vector X. The optimal TSPWT solution
may not exist if the set of time windows cannot be
met due to too short time windows or too long-time
travel distances between customers.
While the variable X(k,1), k
takes the value 1, the corresponding arc between the
nodes (customers) of the graph is part of the
traveller’s cycle. By processing these variables, we
can obtain the order of passing the vertices, i.e.
components of the Matlab vector CYCLE, as the
code shows on lines No. 66 to 79 of the Appendix.
The last flow variables X(k,1), for k
, where , are
specifically the departure times of the vehicle from
the nodes (customers) k at the optimal cycle. The
vector of the departure times is created on line No.
80 of the Appendix.
The real departure time (tStart in Matlab)
of the vehicle from the depot is calculated by a
command on line No. 81. The time (tRet in
Matlab) of the vehicle’s arrival back to the depot
after visiting all customers is calculated by a
command on line No 82. The total duration of the
business trip is calculated on
line No 83.
Table 1. The transformations of variable labels to
Matlab identifiers
Variable
label
Matlab
n
A
b
Aeq
beq
C
D
Variable
label
Matlab
lb
ub
lw
uw
f
kf
m
Variable
label
Matlab
D(i+1,j+1)
X((n+1)*i+j+1,1)
Variable
label
Matlab
C(i+1,j+1)
X(p+1:end,1)
Variable
label
Matlab
X
intcon
t0
X(p+i)
Variable
label
Matlab
tStart
tRet
TotDur
TotDist
For a clear output of the gained solution it is
useful to create vectors of departure and arrival
times (tArr in Matlab) and (tDep in
Matlab) from and to customers with items in the
order of the shortest cycle nodes. It is useful to
extend the vector of departure times by adding
a new first item (the departure time from the
depot). For the last item of the arrival times vector
, it is useful to add (the arrival time to the
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
189
Volume 21, 2024
depot). It is also good to arrange the lower and
upper bounds of time windows
and in the order of the nodes visited. The
vectors , , and are created and
modified using commands on lines No. 84 to 91.
The output time items of all mentioned vectors are
expressed in the clock time format 'hh:mm'. The
total distance (TotDist in Matlab) travelled
during the sales trip is calculated by executing the
commands on lines No. 92 to 99.
5 Application
Our created optimization program was applied to
the delivery of frozen and refrigerated goods from a
central warehouse (vehicle depot) in a company in
the Czech Republic. The optimization of one line
with ten customers is chosen for the presentation.
The GPS coordinates of the vehicle depot are
longitude and latitude .
The customers GPS coordinates, , time window
lower and upper bounds , (clock), and
service times (minutes) are given in Table 2.
Table 2. The customer's GPS coordinates, time
window bounds, and service times
(⁰)
(⁰)
(min)
1
15.7741
49.9841
6:00
8:00
23
2
16.4334
49.9015
10:00
14:00
21
3
14.5875
49.9039
6:00
11:00
22
4
14.8589
49.9939
6:00
8:00
21
5
14.5668
50.0419
6:00
12:00
21
6
16.2156
49.9883
8:00
14:00
22
7
15.3376
49.9463
6:00
8:00
24
8
14.0708
49.9624
6:00
9:00
23
9
14.2490
49.9196
6:00
12:00
21
10
16.2042
49.9052
8:00
13:00
22
The depot and customer locations are plotted in
Figure 1. The data of distances and time distances
bet-ween each two nodes have been obtained
directly from navigation applications. All distances
and time distances are contained in the distance
matrix and time distance matrix in Table 3 and
Table 4.
Fig. 1: Positions of the depot and customers
according to GPS coordinates
Table 3. The distance matrix
Dist.
(m)
j
0
1
2
3
4
5
i
0
0
14172
49463
111603
85370
109170
1
141780
0
59899
102123
75890
99690
2
49454
59734
0
156662
130429
154229
3
112058
102104
156666
0
26711
20196
4
85750
75796
130358
26609
0
24176
5
109477
99523
154085
20138
24130
0
6
27459
40398
22489
137326
111093
134893
7
51454
39296
93858
63847
37614
61414
8
144563
137484
197204
50579
65232
42444
9
135290
129140
183702
35563
53747
33171
10
29914
40194
19780
137122
110889
134689
Dist.
(m)
j
6
7
8
9
10
i
0
27459
50951
143159
137742
29893
1
40571
39392
136367
128162
40329
2
22489
93831
196394
182801
19767
3
137338
63824
50661
34512
137096
4
111030
37516
66234
52748
110788
5
134754
61243
42792
33726
134515
6
0
74495
177058
163465
13220
7
74530
0
103161
89986
74288
8
177876
103944
0
18336
177634
9
164374
90860
18382
0
164132
10
13241
74291
176854
163261
0
Table 4. The time distance matrix
Time
(sec)
j
0
1
2
3
4
5
i
0
0
678
1955
4646
3492
4454
1
678
0
2346
4314
3160
4122
2
1956
2329
0
6382
5227
6189
3
4644
4312
6379
0
1199
945
4
3486
3154
5221
1194
0
1002
5
4446
4114
6180
942
1000
0
6
1238
1643
1050
5695
4541
5503
7
2137
1637
3704
2735
1581
2543
8
6088
5773
7938
2202
2849
1913
9
5642
5325
7391
1470
2211
1467
10
1202
1574
769
5627
4473
5435
Time
(sec)
j
6
7
8
9
10
i
0
1238
2142
6133
5684
1201
1
1662
1637
5810
5352
1593
2
1050
3705
7982
7419
770
3
5695
2733
2208
1488
5625
4
4537
1575
2861
2232
4468
5
5496
2535
1925
1485
5427
6
0
3018
7296
6733
548
7
3020
0
4315
3773
3951
8
7254
4271
0
903
7185
9
6707
3746
906
0
6638
10
548
2950
7227
6664
0
5.1 The Smallest Cycle Duration with
Respecting Time Windows
By running the TSPTW_SOLVER.m function with
the above given input parameters (Table 2 and Table
4), the optimal solution was found. The smallest
duration cycle concerning time windows is given by
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
190
Volume 21, 2024
the sequence of nodes 0-1-7-4-8-9-3-5-10-2-6-0.
The found cycle is drawn in Figure 2. The total time
spent on the delivery route is 8 h 49 min 19 sec and
the corresponding travel distance is 435.650 km.
Fig. 2: The smallest duration cycle of the seller’s
journey with respecting the time windows
Table 5. The arrival and departure times of the
smallest duration cycle concerning time windows
Time
(clock)
Start
The arrived and dep. times of cycle
0-1-7-4-8-9-3-5-10-2-6-0
Re-
turn
i
0
1
7
4
8
9
-
-
6:00
6:00
6:00
6:00
6:00
-
-
6:00
6:50
7:40
8:49
9:27
-
5:48
6:23
7:14
8:01
9:12
9:48
-
-
8:00
8:00
8:00
9:00
12:00
-
i
-
3
5
10
2
6
0
-
6:00
6:00
8:00
10:00
8:00
-
-
10:12
10:50
12:42
13:16
13:55
14:38
-
10:34
11:11
13:04
13:37
14:17
-
-
11:00
12:00
13:00
14:00
15:00
-
The calculated vehicle departure and arrival
times , , and the lower and upper bounds
, of time windows are in Table 5. The time
values shown in the table are rounded to the nearest
minute.
5.2 The Smallest Duration Cycle Without
Time Windows
If all customers do not require any time window,
then there exists n! different cycles that the seller
can use to visit all customers. The mean value of the
travel time after all possible cycles (without service
times) is
. We can find the cycle of
actual minimum travel time. For this purpose, we
modify the previous mixed-integer linear
programming problem by cancelling constraints (5),
(6), and changing (11) to a constraint:
.
(15)
We assume that all customers accept the arrival
of the vehicle at the earliest a. m. By running
this modified M-function, the solution for the cycle
of minimal travel time can be found. This is the
node sequence 0-6-2-10-1-7-4-5-8-9-3-0, and the
corresponding cycle is shown in Figure 3. To solve
the described TSP, the Matlab code published in [4]
can also be used.
Fig. 3: The smallest duration cycle of the
seller’s journey without respecting the time
windows.
The calculated vehicle departure times and
arrival times for each node are included in
Table 6. The time values shown in the table are
rounded to the nearest minute. This business trip
lasts exactly 8 h 36 min 34 sec and the
corresponding travel distance is 419.740 km. This
means in our case that, when respecting the time
windows, the travel time is 12 min and 45 seconds
longer and the distance traveled is 25.910 km
longer.
Table 6. The arrival and departure times of the
smallest duration cycle without time windows
Time
(clock)
Start
The arrived and dep. times of cycle
0-6-2-10-1-7-4-5-8-9-3-0
Re-
turn
i
0
6
2
10
1
7
-
-
6:00
6:39
7:13
8:01
8:51
-
5:39
6:22
7:00
7:35
8:24
9:15
-
i
-
4
5
8
9
3
0
-
9:42
10:19
11:12
11:51
12:36
14:15
-
10:19
10:40
11:35
12:12
12:58
-
5.3 The Smallest Length Cycle without Time
Windows
As is clear from the elements of the used matrices
and , the time distances between the same pairs of
nodes are not exactly proportional to the length
distances. The reason may be, for example, the
different quality of the roads used. Driving time is
shorter on higher quality roads than when using
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
191
Volume 21, 2024
lower class ones. The consequence may be that the
shortest cycle in time may not be the shortest in
length and vice versa.
We can present it in our case. If we optimize our
TSP problem without time windows based on the
use of the distance matrix , we obtain the optimal
solution, given the node sequence 0-5-8-9-3-4-7-1-
10-2-6-0, which is shown in Figure 4. For
optimization we used the code published in [5].
Fig. 4: The minimal length cycle without respecting
time windows
The length of this cycle is 419.441 km and the
total time spent on the delivery route is 8 h 36 min
53 sec. The cycle length is 299 m less than the cycle
length in the previous case, but the travel time is 19
s longer. The last obtained cycle is minimal in terms
of travel distance but not minimal in terms of time.
The calculated vehicle departure times and
arrival times (in hh:mm) for each node are
included in Table 7.
Table 7. The arrival and departure times of the
smallest duration cycle without time windows
Time
(clock)
Start
The arrived and dep. times of cycle
0-5-8-9-3-4-7-1-10-2-6-0
Re-
turn
I
0
5
8
9
3
4
-
-
6:00
6:53
7:31
8:16
8:58
-
:45
6:21
7:16
7:52
8:38
9:19
-
i
-
7
1
10
2
6
0
-
9:45
10:37
11:26
12:01
12:40
13:20
-
10:09
11:00
11:48
12:22
13:02
-
6 Conclusion
This paper proposes a new practical algorithm for
the TSPTW solution for any number of customers,
by using a developed Matlab code. The TSPTW is
formulated as a mixed-integer linear programming
problem with a new approach, which respects the
given matrix of distances, time windows, service
duration times of customers, and constant speed of
the vehicle. The solution lies in minimizing the
vehicle trip duration that services all nodes
(customers). In practice, there is not strict
proportionality between the time and the length of
the path connecting the nodes. Therefore, the
shortest path in time may not be the shortest in
length. The designed objective function using
members with appropriate multiples of departure
times allows minimization of customer departure
times and avoids unnecessary waiting times.
The main output of this paper is a computational
algorithm implemented in an M-function created in
Matlab code, which allows a general solution of the
newly formulated TSPTW problem for any number
n of customers. The created M-function
TSPTW_SOLVER.m is given in the Appendix at
the end of this paper.
The application in section 5 shows one practical
TSPWT solution for 10 customers on an ordinary
PC, this calculation takes 7.23 sec. For comparison,
the optimal solutions for the shortest time and the
shortest distance when no time windows are applied,
are shown.
The M-function TSPTW_SOLVER.m is
practically usable on a regular personal computer
for up to 20 customers. For less than 13 customers,
the calculation takes a few seconds. With a larger
number of customers, the calculation can take tens
of minutes or more. The program was successfully
tested for up to 60 customers.
Entering too short and inappropriately combined
time windows can easily make the problem
unsolvable. It is therefore important to negotiate
with customers and use time windows only when
necessary.
Testing of the created program showed that it is
highly usable in practice. Respecting the results
obtained during the delivery of goods ultimately
enables fuel savings.
The author intends to extend his future research
to solving the vehicle routing problem with time
windows (VRPTW) and respecting the
predetermined priority of customers. In the next
stage of the TSPTW and VRPTW solution, the use
of artificial intelligence is expected.
Acknowledgement:
The submitted manuscript contains a specific model
created by the author for solving TSPTW as an
optimization problem of linear programming (4)
with a specifically designed objective function (13)
and constraints (5) to (12). Service times and time
windows are respected.
The main and unique result of the article is the
transformation of the created model into the Matlab
environment, and the creation of the corresponding
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
192
Volume 21, 2024
M-function, which is applicable for solving TSPTW
for any number of customers. The created M-
function TSPTW_SOLVER.m is in the Appendix.
References:
[1] Bentley, J. J.: Fast Algorithms for Geometric
Traveling Salesman Problems. ORSA Journal
on Computing. Vol. 4, No. 4, 1992, pp. 387–
411.
[2] Davendra, D., Traveling Salesman Problem,
Theory and Applications. InTech, Jerza
Trdina 9, Rijeka, Croatia, 2010.
[3] Gutin, G., and Punnen, A. P., The Traveling
Salesman Problem and Its Variations. New
York: Springer Science+Business Media,
LLC, 2007.
[4] Warren, R. H.: Generating Cyclic Permutati-
ons: Insights to the Traveling Salesman Pro-
blem. WSEAS Transactions on Information
Science and Applications, Vol. 18, 2021, pp.
68–72,
https://doi.org/10.37394/23209.2021.18.9.
[5] Zahradka, J.: The Traveling Salesman
Problem Solution by Mixed Integer Linear
Programming in Matlab Code. Journal of
Applied and Computational Sciences, Vol. 1,
No. 1, 2022, 45–
52.
https://doi.org/10.5281/zenodo.6880928.
[6] Bard, J. F., Kontoravdis, G., and Yu. G.: A
Branch-and-Cut Procedure for the Vehicle
Routing Problem with Time Windows.
Transportation Science, Vol. 36, No. 2, 2002,
pp. 250–269.
[7] Gradle, K. P., and Muley, Y. M.: Travelling
Salesman Problem with MATLAB Program-
ming. International Journal of Advances in
Ap-plied Mathematics and Mechanics. Vol. 2,
No. 3, 2015, pp. 258–266.
[8] Ohlmann, J. W., J. W., and Thomas, B. W.: A
Compressed-Annealing Heuristic for the
Traveling Sales-man Problem with Time
Windows. Informs Journal on Computing,
Vol. 19, No. 1, 2007, pp. 80–90.
https://doi.org/10.1287/ijoc.1050.0145.
[9] Papalitsas, C., Andronikos, T., Ciannakis, K.,
Theocharopoulou, G., and Fanarioti, S.:
A QU-BO Model for the Traveling Salesman
Problem with Time Windows. Algorithms,
Vol. 12, No. 11, 2019, pp.
224.
https://doi.org/10.3390/a12110224.
[10] Sushma, M., Reddy, V.: Finding an Optimal
Path with Hospital Information System Using
GIS-based Network Analysis WSEAS Tran-
sactions on Information Science and Appli-
cations, Vol. 18, 2021, pp. 1–6.
https://doi.org/10.37394/23209.2021.18.1.
[11] MathWorks. Inc., Optimization Toolbox,
User’s Guide. The MathWorks, Inc. Natic,
Massachusetts, United States, 2024,
[Online].
https://uk.mathworks.com/help/optim/index.ht
ml?s_tid=CRUX_lftnav (Accessed Date:
Marcch 11, 2024).
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The sole author of this scientific article
independently conducted and prepared the entire
work from the formulation of the problem to the
final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
This work was funded by the private sources of the
author.
Conflict of Interest
The author has no conflict of interest to declare.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
193
Volume 21, 2024
APPENDIX
TSPTW_SOLVER.m
1: function[X,CYCLE,AllDist,AWT,tStart,tArr]=…
TSPTW_SOLVER(n,D,C,lw,uw,m)
2: Ctemp=C'; Ctemp=Ctemp(:); p=(n+1)*(n+1);
3: for i=1:p
4: if Ctemp(i)==0
5: Ctemp(i)=Inf;
6: end
7: end
8: kf=min(Ctemp)/86400/n; Ctemp=[];
9: CT=C'; f=CT(:); f=[f;kf*ones(n,1)];
10: t0=min(lw-C(1,2:end));
11: A=zeros(n^2,(n+1)^2+n); k=0;
12: for i=1:n
13: for j=1:n
14: if i~=j
15: k=k+1; A(k,(n+1)*i+1+j)=C(i+1,j+1)…
+m(i)+uw(i)-lw(j);
16: A(k,p+i)=1; A(k,p+j)=-1;
17: b(k,1)=uw(i)-lw(j)+m(i) -m(j);
18: end
19: end
20: end
21: for j=1:n
22: k=k+1; A(k,(n+1)*i+1+j)=C(1,j+1)+t0-lw(j);
23: A(k,p+j)=-1; b(k,1)=-lw(j)-m(j);
24: end
25: Aeq=zeros(3*n+3,p+n);
26: for i=1:n+1
27: for j=1:n+1
28: Aeq(i,(i-1)*(n+1)+j)=1;
29: end
30: Aeq(i,(i-1)*(n+1)+i)=0;
31: beq(i,1)=1;
32: end
33: for i=1:n+1
34: for j=1:n+1
35: Aeq(n+1+i,(j-1)*(n+1)+i)=1;
36: end
37: Aeq(n+1+i,(i-1)*(n+1)+i)=0;
38: beq(n+1+i,1)=1;
39: end
40: for i=1:n+1
41: Aeq(2*n+2+i,(i-1)*(n+1)+i)=1;
42: beq(2*n+2+i,1)=0;
43: end
44: lb=zeros(p,1);
45: for i=1:n
46: if lw(i)<t0+C(1,1+i)
47: lb(p+i,1)=t0+C(1,1+i)+m(i);
48: else
49: lb(p+i,1)=lw(i)+m(i);
50: end
51: end
52: k=0;
53: for i=1:n+1
54: for j=1:n+1
55: k=k+1;
56: if i==j
57: ub(k,1)=0;
58: else
59: ub(k,1)=1;
60: end
61: end
62: end
63: ub(p+1:p+n,1)=uw+m;
64: intcon = 1:p; options = optimoptions'intlin…
prog','MaxTime',420,'MaxNodes',3000000);
65: X=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,[],…
options); X=round(X);
66: for i=2:n+1
67: if X(i)==1
68: CYCLE=0; Nok=2; CYCLE(Nok)=i-1;
69: TEST=i; break;
70: end
71: end
72: while TEST~=1
73: for j=1:n+1
74: if X((CYCLE(Nok))*(n+1)+j)==1
75: Nok=Nok+1; CYCLE(Nok)=j-1;
76: TEST=j; break;
77: end
78: end
79: end
80: t=X(p+1:end);
81: tStart=(t(CYCLE(2))-m(CYCLE(2))-C(1,…
CYCLE(2)+1))
82: tRet=t(CYCLE(end-1))+C(CYCLE(end-1)+1,1);
83: TotDur=tRet-tStart;
84: tDep=[tStart,t(CYCLE(2:end-1))];
85: tDep=hours(tDep); tDep.Format='hh:mm';
86: tArr=[t(CYCLE(2:end-1))-(m(CYCLE…
(2:end-1))),tRet];
87: tArr=hours(tArr); tArr.Format='hh:mm';
88: uw=hours(uw(CYCLE(2:end-1)));
89: uw.Format='hh:mm';
90: lw=hours(lw(CYCLE(2:end-1)));
91: lw.Format='hh:mm';
92: TotDist=0;
93: for i=1:(n+1)
94: for j=1:(n+1)
95: if X((n+1)*(i-1)+j)==1
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
194
Volume 21, 2024
96: TotDist=TotDist+D(i,j); break;
97: end
98: end
99: end
WSEAS TRANSACTIONS on INFORMATION SCIENCE and APPLICATIONS
DOI: 10.37394/23209.2024.21.18
Jaromír Zahrádka
E-ISSN: 2224-3402
195
Volume 21, 2024
