@Carl Love   

I think that the use of functions  ln  and  evalf  is impractical,  because it slows down the calculation. As discrepancy for a given partition I take the simple sum of absolute differences.

Compare two approaches for the small set {$ 1 .. 15}:

restart;

ts := time():

S := {$ 1 .. 15}:
AllP := [seq(P, P = Iterator:-SetPartitions(S, [[5, 3]], compile = false))]:
lnp := evalf(ln((`*`(S[]))^(1/3))):
Var := proc (P::(list(set)))
local r; r := evalf(`+`(map(proc (b) options operator, arrow; abs(ln(`*`(b[]))-lnp) end proc, P)[])) end proc: 
sort(AllP, (x, y) ->Var(x) < Var(y) )[1];

time()-ts;

                       [{1, 8, 9, 10, 15}, {2, 5, 7, 12, 13}, {3, 4, 6, 11, 14}]
                                                             322.328

restart;

ts := time():

S := {$ 1 .. 15}: Var := infinity:

T := combinat[choose](S, 5):

for t in T do

S1 := `minus`(S, t):

U := combinat[choose](S1, 5):

for u in U do

S2 := `minus`(S1, u):

Var1 := abs(convert(t, `*`)-convert(u, `*`))+abs(convert(u, `*`)-convert(S2, `*`))+abs(convert(t, `*`)-convert(S2, `*`)):

if Var1 < Var then Var := Var1:

L := {Var, {S2, t, u}} end if

end do end do:

L;

time()-ts;

                                {576, {{1, 4, 13, 14, 15}, {2, 5, 9, 10, 12}, {3, 6, 7, 8, 11}}}
                                                                                 14.516

The results are different, but the discrepancy is actually one and the same, that just shows the non-uniqueness of solution:

[{1, 8, 9, 10, 15}, {2, 5, 7, 12, 13}, {3, 4, 6, 11, 14}];

abs(convert(%[1], `*`)-convert(%[2], `*`))+abs(convert(%[2], `*`)-convert(%[3], `*`))+abs(convert(%[3], `*`)-convert(%[1], `*`));

                                                                576

@Markiyan Hirnyk  In my example, do the rotation by the angle  arcsin (a/2)  about the axis Oy and the homothety with coefficient  sqrt (3) and the center at the origin. Then you will get the same coordinates as in Wiki.

@Markiyan Hirnyk 

1. Both procedures are very simple conceptually and based on brute force technique of two or three vertices and then checking (using geometric properties listed above). So I do not think there is a need for separate comments.

2. To be more precise, that used O (n ^ 3) algorithm, because  binomial(n, 3)=n*(n-1)*(n-2)/6 . Of course, in terms of time, these procedures are not optimal and appropriate for relatively small polygons (n<100). But their advantages in simplicity and exact solutions. Therefore they can be used in the proofs in geometry problems.

3. I do not agree with this statement "Exact solutions and close-form-expressions are from somewhat outdated math" . The advantage of an exact solution at least in the fact that from it is easy to get  an approximate solution with any precision, but not vice versa.

@toandhsp  I was looking for it of right triangles with legs   x+3  and  y+3  and the hypotenuse  x+y . Radius of the inscribed circle is  3  for any such triangle.

isolve((x+3)^2+(y+3)^2=(x+y)^2):
op(select(s->rhs(s[1])>0 and rhs(s[2])>0,[%]));

{x = 4, y = 21}, {x = 5, y = 12}, {x = 6, y = 9}, {x = 9, y = 6}, {x = 12, y = 5}, {x = 21, y = 4}

@toandhsp   [[4,-2], [4,7], [-8,-2]]

see  http://www.mapleprimes.com/posts/200122-The-Largest-Incircle-And-The-Smallest-Covering-Circle-

@Markiyan Hirnyk  Nice and simple solution. Now I am preparing two procedures in which  both problems will be solved exactly (symbolically). For example, the exact answer to this example is

[[57/22, 43/22], 1/242*1625^(1/2)*242^(1/2)];
evalf[15](%);

@sainimu78  Syntax error! Under the symbol  ∑  the command  add  should not be used.

@isbah It is impossible to obtain solution explicitly in the form of formulas. The reason is that the solutions are the roots of the polynomial of high degree  RootOf(...) . It is well known that if the degree of a polynomial more than 4, then the roots are not expressed as a function of its coefficients. But the procedure

f:=a->fsolve({0.5704+0.5211*exp(-0.6569*x)-0.03849*y=a, 0.561+0.5284*exp(-0.6479*x)-0.03929*y=a}):

can be used also successfully as  explicit formulas, for examples, to find the values ​​at certain points, build plots, etc.

@Carl Love  In initial task, we don't need use these commands, because they slow down the program and we can do the work without them.

Example:

A:=[1,4,8]: B:=[8,9,12]:
type(sqrt(A[1]^2+A[2]^2+A[3]^2), posint);
type(sqrt(B[1]^2+B[2]^2+B[3]^2), posint);
type(sqrt((A[1]-B[1])^2+(A[2]-B[2])^2+(A[3]-B[3])^2), posint);

                                               true

                                               true

                                               false

The main disadvantages of the source code - it is slow and duplicates occurs.

It is too difficult and it is not well suited for this forum. Look at the links found by Google

https://www.google.ru/#newwindow=1&q=Maximum+area+rectangle+inscribed+in+a+convex+polygon

@Markiyan Hirnyk  Of course my reasoning is incorrect. Thank you for this error message.

@Markiyan Hirnyk  I meant that, for example, both numbers  F1(x0, y0, z0))  and  F1(x5, y5, z5)  are positive or both ones are negative numbers. In both cases should be  F1(x0, y0, z0))*F1(x5, y5, z5)>0 . Function  sign  we have not use.

Unfortunately, Maple does not solve even the simplest equations containing the function  floor :

solve(floor(x) = 1, x);

                RootOf(floor(_Z)-1)

@Markiyan Hirnyk  If you remove the empty line   [``,``,``,``,``]  then there are no problems.