The reason why is explained in the attached file.
(Note that the main problem comes from your defining four constraints whereas there can't be more than two -as  you write-,  see the attached file).

Parameter whose numeric values have to be given are [Pu, U, a, beta, eta, tau0, upsilon, varphi, w].

Once done you will get 3 solutions which depend on {Cm, Cr, lambda, varepsilon}.  
The fact the inside point is a maximum or a minimum depends on lambda  and varepsilon and it is impossible to conclude as long as (1) they are not given numeric values or (2) some assumption on lambda (see relation (16) ).

Q_P1_Optimum_condition(1)_mmcdara.mw

Watchout: the last results (lambda-varepsilon relation  in equation (16) ) and plots:-inequal depend on the seed used and will change if you execute the data := ...; command without restarting the code (for the seed changes).

Indeed, isolate help page says "If the first argument is not an equation, then the equation eqn = 0; is assumed."

Here are very simple examples where isolate does nor return the (likely) expected result

isolate(p+b <= d, p) 
                          p <= -b + d
isolate(a*p+b <= d, p); # the solution depends on the sign of a
                         a p <= -b + d
# if a > 0
assume(a > 0):
isolate(a*p+b <= d, p);

                               -b + d
                          p <= ------
                                 a   
# if a < 0... attempt 1
a := 'a': 
assume(a < 0): 
isolate(a*p+b <= d, p);     
                    p a <= -b + d
# if a < 0... attempt 2
a := 'a': b := 'b': d := 'd':
assume(a < 0, b > 0, d > b):
isolate(a*p+b <= d, p);
                          p a <= d - b
# and even simpler
​​​​​​​isolate(-3*p <= 2, p);
                           -3 p <= 2

So I advice you using solve instead
 

Download Q_solve.mw

D is a protected name which stands for Derivative, search for D in the help page.

Insert the line  local D; before line # Define coefficients for RSM (adjust as needed)
RSM_model.mw

I am right with @Rouben Rostamian : the area of degenarate ("flat") triangle is 0 and the collinearity condition is not mandatory tp compute the area.

Nevertheless some other constructions based on the triangle (excircle, circumcircle, ...) require that the triangle is not degenerate. I fell like the non collinearity condition is introduced from the start (the triangle command) to avoid asking him each time it becomes mandatory.

The attached file provides a way to assume the non collinearity condition without copying and pasting it manually as you did.
Triangle_given_condition.mw

You wrote "Under what conditions on the parameters does the lower bound condition become violated while the upper bound condition is also violated?"

If you maximize/minimize some function under inequality constraints defining a multidimensional domain D, and if some of them arestrict inequalities, either the extremum lies in the open D, either it belongs to the boundary of D. But the constraints are never violated.

Whatever, does this suit you?
Q_Const_mmcdara.mw

By the way, you use both U and U[0], which is dangerous and may lead to wrong results.
So I changed U[0] into U__0.

The 'optimize' option in CodeGeneration doesn't seem capable to get the representation you desire (with Maple 2015.2)

So I propose you a code (not written as a procedure in order you can see easili how it works without using the debugger) which could suit you.

Basically it "extracts" all the indeterminates of any type from an expression (possibly a matrix of expressions), iteratively build rewritting rules (your t1 = a*b*c or t2 = d+e) and next apply them to the expression (M) to "reduce" it.
There are several possibilities, which depend on the order in which you consider those indeterminates (a point that I did not realy investigate) and on how much rewritting rules you use.
For instance if M is a 2-by-2 matrix whose elements are very lengthy expressions, then a hypothetical CodeGeneration output like this one would give the most compact form of some "optimized M"

t1 = M[1, 1]
t2 = M[1, 2]
t3 = M[2, 1]
t4 = M[2, 2]
M[1, 1] = t1
M[1, 2] = t2
M[2, 1] = t3
M[2, 2] = t4

But this to the detriment of lenghty expressions of t1, ..., t4.

So I defined a complexity measure to rank those possibilities according to the "complexity" of the transformed expression of M.

As you will see adding more rules does not necessarily over simplify theis expression.
For instance, according to the complexity measure I choosed, here is a simpler representation of M:

t1 = d+e
t2 = a*b*c
t3 = t1*t2
M[1] = t3
M[2] = t2-t1

A more complicated example than yours is also presented and a new complexity measure, which accounts both for the rewritting rules and the final representation of M is proposed.
Optimized_code_0.mw

NOTE: this code will likely fail if tested on more complex situations.

that could interest you:

1. How to define a truncated normal random variable and sample it efficiently?
   (drawing a sample of size 40 with the acceptation/rejection code you write does is rather fast... but what if you want to draw a sample of size 10⁴, or more?
    
2. Tally: function to get frequencies (a TallyInto also exists)
    
3. FrequencyTable: a function to get more informations about frequencies (this point is related to section Ajouter les fréquences & fréquences cumulées).
    
4. Histogram(..., discrete=true) to get a bargraph.
    
5. ... and maybe a few others
    
6. About your first "tableau"  (Mettre les données brutes dans un tableau), you could use DocumentTools:-Layout (see the corresponding help page ant the exemples provided) to ebhance ir. For ewample by giving rows names which refer to the distributions those numbers are drawn from.
   More of this you could make nu and varphi appear as greek letters directly.
    

Truncated_Normal_Distribution.mw

You will find in worksheet Generating_random_number_efficiently.mw that the maon cause of the inefficiency of your code is that you generate each random number once at the time.
I show how to generate all of them in a row (even for a truncated law) and thus enhance dramatically the performance of the code.

The first 52 roots with |error| < 1e-7: Find_roots.mw

Apply yourself the same stuff for case 2.

Note: the index in your for loop should not be varepsilon, that is too dangerous: look to this

restart:
f := a*varepsilon^2:

for varepsilon from 1 to 2 do
  f;
end do;

f;  #the value of the index being equal to 3
                               a
                              4 a
                              9 a
#-----------------------------------------------------
restart:
f := a*varepsilon^2:

for val from 1 to 2 do
  eval(f, varepsilon=val)
end do;

f;
                               a
                              4 a
                                     2
                         a varepsilon 

See Question_New_mmcdara.mw for case 1... the plot you aske for will not give you any valuable informations given the ranges of variations of Pc, p1 and p2.

More of this analytic solutions are easily got:

1. Your first maximization problem can be solved formally using lagragian multipliers.
2. Your second maximization problem is even simpler as it reduces to maximize a function (TP) which is quadratic wrt p1 and p2.
   Taking its first partial derivatives show that diff(TP, p1) = A(varepsilon)*p1+B(varepsilon) ... but diff(TP, p2) = C*p2+D: so the "optimum" p2 does not depend on varepsilon.
   Plotting "optimum" p2 against varepsilon is thus useless.
   
    

   >  >  >  >  >  >  >  >  >  >  >  >  >  >  >  ANALYTIC SOLUTION >  obj   := TRC: cstr  := [C_1[], C_2[]]: H     := lambda__1*(lhs(cstr[1])-Pc), lambda__2*(Pc-rhs(cstr[2])): lag   := obj + add(H): sols  := [solve([diff(lag, Pc), H], {Pc, lambda__1, lambda__2})]: sols  := map(s -> if eval(lambda__1^2+lambda__2^2, s) = 0 then NULL else s end if, sols): objs  := map(s -> eval(diff(obj, Pc), s), sols): pairs := zip((s, t) -> [eval(Pc, s), t], sols, objs); print():print(): TP_1   := eval(L_wout(p1, p2), sols[1]): TP     := collect(subs(DATA, TP_1), [p1, p2]): sol    := solve(diff~(TP, [p1, p2]), {p1, p2}): triple := eval([p1, p2, TP], sol); (1) >  with(plots): >  display(   plot([pairs[1][], varepsilon=0..2], color=blue, linestyle=3, labels=[Pc, Max(TR)])   ,   seq(     textplot([eval(pairs[1], varepsilon=val)[], typeset(varepsilon=nprintf("%1.1f", val))], align={right})     , val=0..2, 0.5   )   ,   seq(     pointplot([eval(pairs[1], varepsilon=val)], symbol=solidcircle, symbolsize=10, color=blue)     , val=0..2, 0.5   )   , view=[1400..1700, 10^5..1.6*10^5]   , size=[500, 400] ); display(   plot([triple[[1, 3]][], varepsilon=0..2], color=blue, linestyle=3, labels=[p1, Max('TP')])   ,   seq(     textplot([eval(triple[[1, 3]], varepsilon=val)[], typeset(varepsilon=nprintf("%1.1f", val))], align={right})     , val=0..2, 0.5   )   ,   seq(     pointplot([eval(triple[[1, 3]], varepsilon=val)], symbol=solidcircle, symbolsize=10, color=blue)     , val=0..2, 0.5   )   , size=[500, 400] )   >  display(   plot([pairs[1][], varepsilon=0..2], color=blue, linestyle=3, labels=[Pc, Max(TR)])   ,   plot([pairs[2][], varepsilon=0..2], color=red, linestyle=3)   ,   seq(     textplot([eval(pairs[1], varepsilon=val)[], nprintf("%1.1f", val)], align={right})     , val=0..2, 0.5   )   ,   seq(     pointplot([eval(pairs[1], varepsilon=val)], symbol=solidcircle, symbolsize=10, color=blue)     , val=0..2, 0.5   )   , size=[500, 400] ); >  display(   plot([pairs[1][1], triple[1], varepsilon=0..2], color=blue, linestyle=3, labels=[Pc, p1])   ,   seq(     textplot([eval([pairs[1][1], triple[1]], varepsilon=val)[], nprintf("%1.1f", val)], align={right})     , val=0..2, 0.5   )   ,   seq(     pointplot([eval([pairs[1][1], triple[1]], varepsilon=val)], symbol=solidcircle, symbolsize=10, color=blue)     , val=0..2, 0.5   )   , view=[1400..1700, 19000..25000]   , size=[500, 400] ); >  dualaxisplot(   plot(     pairs[1][1]     , varepsilon=0..2     , color=red     , legend=Pc     , labels=[typeset(varepsilon), Pc]   )   ,   plot(     triple[1]     , varepsilon=0..2     , color=blue     , legend=p1     , labels=[typeset(varepsilon), p1]   )   , size=[500, 400] ); >

    

   Question_New_mmcdara(analytic).mw

I do not pretend the way I take i the attached file is perfect, this is just an idea based on the use of a filter to remove the the discontinuities while damping them.
So the plot is obviously wrong close to the discontinuities, but no longer contains those vertical rays which seem to bother you.
plots-long_term_mmcdara.mw

Note: the last plot is very far from being perfect and a lot of work should be necessary to display what I caled the "discontinuity wall" in a smooth way.

@GFY 

kk(3, m) posseses 4 branches and kk(3, 200)[1] returns a value which is not the one plot(kk(3,m)[1], m=200..600) displays at m=200. Right?
Simply because plot does not select the same branch for some reason I don't know.
So this is a workaround:

 

Download gra423_Omega_Maple2015.mw

Simply write Sum instead of sum (Maple 2015)

L := limit((Sum(Sum(j/(j^2+k^2), k = 1 .. n), j = 1 .. n))/n, n = infinity)
                         1         1   
                         - ln(2) + - Pi
                         2         4   

Wait for a detailed explanation by some Maple expert (@vv maybe?)
 

Explanation for the first plot alone (I guess all the other errors are of the same nature)

Download Understand_your_errors.mw

And now two examples of how you could proceed
 

Download What_you_could_do.mw

Insertion can be made using ImageTools: insert.mw

@JaneCherrytree 

Do you use the name PDF on purpose?
I mean, does the variable PDF represent the Probability Density Function of some random variable?

If it is so PDF must be nonnegative over (x, y) in (0..1)^2... which is not be the case given the ranges of the five parameters VeX, VeY, VeZ, alphaX, alphaY that you provide ^(#).

Next, there is absolutely no chance that a formal expression of the (CDF?) does exist, restricting you to numerical constructions alone.

See the attached file IntNum.mw

If this suits you feel free to ask some extra developments

^(#) Add something like this at the end of the attached file to see that the proportion of instances which leed to a negative scaling (meaning PDF is not a Probability Density Function) is about 50%

NegativeScalings := 0:
for rep from 1 to 100 do
  data := [seq(rand(-1. .. 1.)(), i=1..3), seq(rand(0. .. 1.)(), i=1..2)]:
  instances := params =~ data;
  scaling := IntNum(instances, 0..1, 0..1);
  if scaling < 0 then NegativeScalings := NegativeScalings + 1 end if:
end do:

NegativeScalings;
                               54