I took the information from the StackExchange article that you referenced and use it to write the procedure LinearizeAbsObj below. Then I used that to resolve the same problem as an ILP with Optimization:-LPSolve. Now the solution is almost instantaneous (16 ms).
 

Download Taxi_metric_ILP.mw

If A is any real square nonsingular matrix, then A times A-transpose will be symmetric positive definite (i.e., all eigenvalues are positive), which is suitable for Cholesky. The probability that a random matrix is singular is infintesimal, so just take any random real square matrix and muliply it by its transpose. In Maple, the short form for A-transpose is A^+ or A^%T (the former syntax is newer).

A:= LinearAlgebra:-RandomMatrix(6$2):
M:= A.A^+:
L:= LinearAlgebra:-LUDecomposition(M, method= Cholesky);

The issue isn't so much that A is nonsquare; it's that any Matrix/Vector/Array operation (not just LinearSolve) that uses option inplace is not allowed to change the size of the container which is supposedly being operated on "in place". So, it fills unused spots with 0. I suppose that a more-logical option might be to use NULL as the fill value, but that wouldn't work if the container were declared sparse.

Addressing your second question: There any many possibilities. Are you looking for an exact solution? Is the system overdetermined (as one would ordinarily expect when rows > columns)? If so, are you looking for a solution in the least-squares sense?

Using matrices/vectors to store the variables Y[i,j] and X[j] makes it easier to enter the objective and the constraints and to display the final results.
 

Download logistics_ILP.mw

To see the prime factorization in its raw list form, use ifactors instead of ifactor:

ifactors(8650368);

     [1, [[2, 7], [3, 3], [2503, 1]]]

ifactors(8650368)[2];

     [[2, 7], [3, 3], [2503, 1]]

To see that expressed as a product of powers, do

InertForm:-Typeset(`%*`((p-> ``(p[1])^p[2])~(ifactors(8650368)[2])[]));

     (2)^7 * (3)^3 * (2503)

This isn't really a question about Maple because any operating system should have utilities that track and log the usage of any and every program. Any half-way competent sysadmin knows how to access this information . 

Assuming that you're solving an integer linear program (ILP) via Optimization:-LPSolve, you would specify that the variables x[2,3] and x[1,4] (for example) can only take the values 0 or  1 by including the option binaryvariables= {x[2,3], x[1,4]} in the command. See ?LPSolve.

Every time that you have L by itself on the left side of := it changes what L is. So, after your second line, L is no longer an Array, it's a sequence. Here's how to add a new element position:

L:= Array([one, two, three, four, five]):
L(numelems(L)+1):= six:

The spelled out numbers could be anything, and the left side of the second line could be simply L(6).

Are you sure that the original had a z-axis label? A basic plot3d command such as

plot3d(x^2+y^2, x= -2..2, y= -2..2);

doesn't produce a z-axis label unless you explicitly include one via, for example,

plot3d(x^2+y^2, x= -2..2, y= -2..2, labels= ['x', 'y', 'z']);

If that doesn't cure your problem, please upload the code that generates the plot.

fsolve returns a sequence of roots (when solving a polynomial) rather than a set of roots, but select expects a set. So, replace fsolve(r, x[2]) with {fsolve(r, x[2])}. When you do that, you're just going to get the empty set on each loop iteration because r has no roots between 0 and 1; indeed, it has no real roots at all. To see the complex roots, use fsolve(r, x[2], complex) without the select.

The code in the attached file runs without error, so I don't know what to correct. Perhaps you simply mean that you want to see the results? For that use

eval(Result_Eva_GB_X3_X2);

after the loop.

The results of code inside loops is ordinarily not shown. That can be adjusted by setting the value of printlevel (see ?printlevel). One needs to be careful with this because it's very easy for the amount of information displayed to become overwhelming, and once it has started displaying there's usually no way to stop it other than killing the whole kernel.

solve({coeffs(expand(u-f), x)});

However, you have five variables and only four equations, so you won't get a complete solution.

In the following three function definitions

Ty:= y0-> T*diff((r-> rhs(r[-1]))@pds:-value(T(x, y), x = 0), [1], [y0]); 
Ty2:= y0-> T*diff((r-> rhs(r[-1]))@pds2:-value(T(x, y), x = 0), [1], [y0]); 
Ty3:= y0-> T*diff((r-> rhs(r[-1]))@pds3:-value(T(x, y), x = 0), [1], [y0]);

you need to use numeric differentiation, fdiff, rather than diff (which is only intended for symbolic differentiation). I don't have time to go into the details right now, but it looks like you were using fdiff syntax anyway.

You may also need to reduce the precision of the numeric differentiation for it to return results. Four digits of precision should be enough for any plot. I think that numpoints= 1000 is excessive given the number of mesh points used by pdsolve. You did say spacestep = .25. And I'm not saying that your mesh should be ~1000 points in the space dimension! That would probably take too long to solve.

Also, to ensure that you get multiplication, make sure that you enter "T space bar diff" or T*diff. In the last two of the above three cases, you actually had the single word Tdiff, which is meaningless, and I corrected it. But I'm fairly certain that your attempt to use T as a functional operator like that is incorrect. T and T(x,y) are just symbols; you need to use pds:-value again to access/create a numeric functional operator that corresponds to T.

You can either use

if Equal(Row(A,2), ZeroVector[row](6)) then ....

or (my preference)

if andmap(`=`, A[2, ..], 0) then ....

You can localize the value Digits to a particular call to evalf like this:

evalf[1500](Pi);

This avoids the need to explicitly set the value of Digits.

A lot of your computation can be avoided: If x is a float, then op(1, x) is its base-10 integer mantissa. For example, op(1, 3.14) is 314. Integers can be easily split into their digits with convert(..., base, 10). This works for any base if you just change the 10.

Putting it all together, your computation can be reduced to the single line:

Statistics:-Histogram(convert(op(1, evalf[1500](Pi)), base, 10), discrete, thickness= 9);

The appearance of the plot can be improved with some options:

Statistics:-Histogram(
   convert(op(1, evalf[1500](Pi)), base, 10),
   discrete, thickness= 9, view= [-1..10, 0..0.12], axis[1]= [tickmarks= [$0..9]]
);