I tried to apply your procedures (in Maple 2017.3 and 2016.2) to the same example and got a strange error:

restart;
ranges:=proc(simpledom::list(relation), X::list(name))
local rez:=NULL, r,z,k,r1,r2;
if nops(simpledom)<>2*nops(X) then error "Domain not simple!" fi;
for k to nops(X) do    r1,r2:=simpledom[2*k-1..2*k][]; z:=X[k];
  if   rhs(r1)=z and lhs(r2)=z then rez:=z=lhs(r1)..rhs(r2),rez; #a<z,z<b
  elif lhs(r1)=z and rhs(r2)=z then rez:=z=lhs(r2)..rhs(r1),rez  #z<b,a<z
  else error "Strange order in a simple domain" fi
od;
rez
end proc:

MultiIntPoly:=proc(f, rels::list(relation(ratpoly)), X::list(name))
local r,rez,sol,irr,wirr, rels1, w;
irr:=[indets(rels,{function,realcons^realcons})[]];
wirr:=[seq(w[i],i=1..nops(irr))];
rels1:=eval(rels, irr=~wirr);
sol:=SolveTools:-SemiAlgebraic(rels1,X,parameters=wirr):
sol:=remove(hastype, eval(sol,wirr=~irr), `=`); 
add(Int(f,ranges(r,X)),r=sol)
end proc:

R := [ (x-4)^2+y^2 <= 25, x^2+(y-3)^2 >= 9, (x+sqrt(7))^2+y^2 >= 16 ];
MultiIntPoly(1, R, [x,y]): 
radnormal(simplify(value(%)));

@mapleatha  In this case use  Standard GUI  in which there are no any parentheses around the exponent.

In Classic interface, you can do so:

expa := 2: 
exp(cat(expa,` t`));

Edit.

@vv  For better perception, I did everything on one plot, eliminated some repetitions in the code and added the  fieldstrength = log   option that equalizes the lengths of the arrows:

restart;
with(plots): with(plottools):
poly:=x>=0, y<=1, y>=-1, y<=4-4*x, y>=4*x-4:
F,G:=piecewise~(`and`(poly), [x^2+y^2,3*sin(x)-y^2], undefined)[]:
display(polygon([[0, -1], [3/4, -1], [1, 0], [3/4, 1], [0, 1]], style = line, thickness = 2), fieldplot([F,G], x =0 .. 1, y = -1 .. 1, arrows = SLIM, color = x, fieldstrength = log));

@nunavutpanther If  X  is a list then  add(X)  only works in the latest versions of Maple. Replace   add(X)  with  `+`(X[ ])

@digerdiga 

expression:=(a*b^2*c*d+a^2*b*c*e)*x^2 + (a*b+a^2*e)*x:
map(factor, expression);
 

@Joe Riel   From the condition of the problem  "...carpet  covers exactly 24 plates and at the same time each vertical and each horizontal row of plates in the living room contains 
exactly 4 plates covered with carpet..." it follows that the room size is exactly  6 by 6

In the general case, we also have 5 roots. So what's the problem?

f:=1/4*(-1/4+alpha*(-1+b)*e^2+((1-b)*alpha+b/4)*e)^2/(-1+e)^2/(b*e-1)/alpha/e/(-1+b):
diff(f, e):
solve(%, e):
nops([%]);
                                                           5

@Markiyan Hirnyk  Thank you.

I can not agree with your statement  "This approach is rather a handicraft than a mechanization because a parametrization of the boundary of a domain  is needed."  And if some region is initially defined by the equation of its boundary? In addition, not every region can be written as the intersection of individual regions. For example, a nonconvex polygon can not be defined in this way, because the intersection of half-planes is always convex.

@maxburakho  See Help on  ?fsolve,detail . I did not find there any options that would allow me to force a solution method. In fact, fsolve command itself chooses the most suitable method for calculating the root. For your case, the system is actually solved exactly (not by Newton method), solve command gives the same result. Probably you will have to write code for solving the system yourself, using Newton's method or to look for something similar on the Internet.

@lucaud  I did not understand your request regarding volume. The volume is already calculated in Carl's answer (2 ways) and in my answer (1 way).

When creating complex animations, it is more convenient and quicker to use a procedure to create one frame for an arbitrary value of the animation parameter, and then output successive frames with plots:-animate command.

Here is this technique for your example. I also increased the number of frames for better perception and change the color, otherwise it looks rather gloomy:

restart;
Australopithecus := [[75, 25], [97, 30], [93, 40], [93, 45], [83, 50], [80, 55], [79, 60], [81, 73], [74, 76], [68, 81], [60, 82], [50, 83], [40, 80], [30, 71], [25, 60], [24, 50], [25, 37], [15, 33], [10, 30], [45, 10], [55, 16], [65, 10], [80, 8], [93, 14], [96, 24]]:
man := [[95, 39], [113, 40], [111, 47], [118, 53], [113, 62], [109, 72], [112, 88], [112, 95], [107, 112], [99, 117], [85, 122], [72, 122], [49, 117], [36, 104], [31, 78], [39, 52], [43, 43], [44, 34], [39, 16], [73, 3], [81, 17], [98, 14], [105, 17], [104, 26], [111, 33]]:

morph := proc(poly1, poly2, t)
uses plots, plottools; 
if nops(poly1) <> nops(poly2) then ERROR("mensaje.") end if; 
display(polygon((1-t)*poly1+t*poly2, color = gold), axes = none); 
end proc:

plots:-animate(morph, [Australopithecus, man, t], t = 0 .. 1, frames = 60, paraminfo = false);

Formulate more exactly what you want: volume, picture or something else. Specify these 2 solids here.

Here is a procedure based on vv's idea:

restart;
LR:=proc(X::list)
local X1, y0, X0, L, x0, Y0, y1;
X1:=simplex:-convexhull(X);
y0:=min(seq(x[2], x=X1));
X0:=sort(select(t->t[2]=y0, X1),key = (x->x[1]));
L:=select(t->t[1]>=X0[-1,1], X1);
x0:=max(seq(l[1], l=L));
Y0:=sort(select(t->t[1]=x0, L),key = (x->x[2]));
y1:=min(seq(y[2], y=Y0));
remove(t->t[2]>y1, L);
end proc:

Example of use:
LR([[1,1],[1,2],[1,3],[2,1],[3,2],[3,6],[4,2],[4,6],[5,4],[5,5]]);
                  [[2, 1], [4, 2], [5, 4]]

Visualization of this example:
A:=plot([[1,1],[1,2],[1,3],[3,2],[3,6],[4,6],[5,5]], style=point, color=red, symbol=circle, symbolsize=20):
B:=plot([[2,1],[4,2],[5,4]], style=point, color=red, symbol=solidcircle, symbolsize=20):
P1:=plot([[2,1],[4,2],[5,4]], color=blue, thickness=2):
P2:=plot([[5,4],[5,5],[4,6],[3,6],[1,3],[1,1],[2,1]], linestyle=3, color=blue):
plots:-display(A,B,P1,P2, view=[0..5,0..6], scaling=constrained);

The code was edited.

@Earl  As a point corresponding to the initial value of the parameter t=0, we can take any point of this circle. As you can see from the text of the code, first I found the center of this circle (the list C=[1,1/2,1/2]), so the plane passing through the center intersects the circle at least in two points. Solving the corresponding system, I found 2 points and as a point corresponding to t = 0, I took a point with a larger first coordinate.

@max125  It's like a staircase:

plot3d(trunc(x), x=0.99..4.99,y=0..5);