@vv 

I'm not sure the explanation you proposed in add, floats, and Kahan sum - MaplePrimes applies here.
Indeed, add cannot minimize a numeric error because not all the terms are numeric, more of  this, replacing the lines

Y := G(1, 2, -1)~(X) +  Sample(Normal(0, 0.1), N)^+:

by 

Y := G(1, 2, -1)~(X) +  LinearAlgebra:-RandomVector(N, generator=-10..10) /~ 10:

in order not to have floats but rationals, still shows the same behaviour as with floats.
See sum_changes_ranks_rational.mw.

The strange behaviour still holds in sum_changes_ranks_rational_permutation.mw  (same "noise" added after having been shuffled), and even it does in this extremely simple example sum_changes_ranks_simplest_non_random.mw.

@dharr 

It is a shame that approxsoln requires a guess for all the functions.
For instance the guess G(eta) ~ exp(-eta) seems quite natural and it would be great if approxsoln = [G(eta) = exp(-eta)] was authorized.

I vote up

@KIRAN SAJJAN 

Open a new worksheet, copy-paste this command

?(dsolve/numeric_bvp,advanced)

and execute it.

The error (in dsolve/numeric/bvp) initial Newton iteration is not converging is explained and some cures are proposed... which doe not mean there are easy to apply.
Generally the continuation method is used to ease convergence: it is quite simple to apply when you have a single ode, but you often have to do a lot of trials and fails  when you deal with an ode system (which ode(s) or (and) bc(s) is (are)  responsible of the non convergence).
I spent more than one hour but wasn't capable to fix this non-convergence problem.
Generally a good physical knowledge of the system to solve may help identifying the non-convergencesource (some "boundary layer" somewhere).

Another remedy consists in telling dsolve that you know, or may infer, an approcimate solution.
See here ?(dsolve/numeric_bvp) and search for approxsoln.

I can't do more for you.

@acer 

It would be interesting to know where those a1, ..., a8 expressions come from.
Do they have some "physical" meaning (which I doubt of given the somplexity of some of them) or do they come from a visual observation of ra[1] (which I doubt too as, for instance, a1 doas not appear in ra[1]) ?

At the end, what does the OP really want to achieve?

Can you upload your workshhet (no special characters in its name) using big green arrow in the menubar?

@MAXR 

I thought I explained you with enough details that your design of experiment (a so called composite central design) wasn't consistent with the model to fit.

Your last worksheet proves me I was wrong on all the points.

In my previous comment I also explained that you not having any random error elsewhere the discrepancy between the predictions from the fitted rduced model (quadratic terms removed) and the "true" model (the one which generates the pseudo-observations) was a deterministic function (which 
I denoted H(A², B², C², D²) ).
To beeven clearer what you are doing here has absolutely nothing yo fo with Statistics (randomness is nowhere) but is simply linear algebra.
Indeed the model can be derived out of any probabilistic model of the observation-regressors assumpions. (see my initial answer , equation (10) wher I just minimize thr L2 norm of Nu-G(A, B, C, D) wrt the parameters of G(A, B, C, D)).

Statistics enters the picture only once you have written such probabilistic method:

1. I have quantity Y which depends on some other quantities X through the relation Y = F(X).
2. I observe Y through some device whose, I know, produces a random  observation error  E(X).
3. I assume that the expecation of E(X) is everywhere equal to 0 (just to simplify the sequel).
4. Then the expectation of the observed value of Y at some point X = x is:
                    Expectation(Y_obs | X = x) = F(x) + 0 
5. The model F(X) contains some parameters P (so it's better to write now F(X; P) instead of F(X)) whose values are unknown and that I want to asses given a collection of observations of Y at different locations x₁, ...x_N.
6. To do this I assume that:
   1. The variance of E(X) is finite (mandatory assumption).
   2. The variance of E(X) is independent of X (not a mandatory condition but a common one).
   3. The expectation of E(X =x)_*E(X = x') is equal to 0 for every couple (x, x') (not a mandatory condition but a common one).
7. Given assumptions at points 2, 3, 6 the quantity
              S²(p) = Sum(r_n(p),  i=1..N) / (N-1),  
   where r_n(p) = Y_obs | X = x_n) - F(x_n; P=p) is the nth residu when parameters P have values p, is an estimation of the variance of E(X).
8.  My goal is to determine the value p* of P which minimizes S²(p).
9. In the case where F(X; P) is a linear form wrt P, S²(p) is a quadratic form wrt P and thevalue of p* can be obtained through simple linear algebra.

What you missed in what to do are points 2, 3, 4, 6, 7. So you do not do probabilities nor statistcs and it is pure nonsense to even speak about ANOVA or significativity of the values p*: because you have no randomness all those values are exact and their covariance matrix is exactly equal to 0.

(
Even if Fit or LinearFit say otherwise, but this is become they both assumed you used the probabilistic model defined by points 1 to 6. To be more rough: the assume that you do the things in the state of the art and that you're not doing anything (just as you did
).
 

You want me to help you "to construct the tables in the same format as Tables 6 and 7." ... but did you, a single second, give a look toTable 5?
Had you did so you would have seen the author(s) used a randomized replicated design and that the values of the response at replicates of the same location are not the same, very likely denoting the presence of an observation error somewhere else detailed in the paper (for instance the 5 replicate of the central points give 5 different values of the response).
If you want to "to construct the tables in the same format as Tables 6 and 7." then begin by constructing a table which looks like Table 5.


As you keep using a non observable model despite what I said in my initial answer

(
Observe that in Table 5 the runs 10, 12, 13, 15, 16, 18 are the axial points I told you about which enable assessing the quadratic effects. Thus why the model with quadratic terms the paper presents is identifiable (ie of full ra,k).
Why do you persist in not introducing those axial runs to avoid having a non observable model?
)

read this Final_answer.mw worksheet to see how Fit and LinearFit behave when the model contains non identifiable terms: as LinearFit fires a warning but nevertheless gives a result, Fit fires an error: personally, I prefer to receive an error and no result than a warning and a result whose quality is legitimately questionable.

I won't answer any more of your questions as long as you persist in your positions and refuse to understand your mistakes.

ß@MAXR 

But if you want, as uou say, "assess the significance of model terms in a regression model" why not simply use pvalues?

ANOVA is used when you have, for instance, N observations depending on, for instance, two factors, let's say A and B;
and that you want to assess the effect of A, B, A*B. I don't see the need to use this here.

More of this you must be aware that assessing  the significance of model terms is, in the special case you built, something quite delicate because the underlying assumptions those significance tests require are not met;

• When you do OLS (Ordinary Least Squares) regression you start writting that the expectation of Nu(A, B, C, D) is the model you fit and that the difference between the realisation of the fitted model and the observations (vector Nu) comes from the presence of an additive error E(A, B, C, D).
• This error E(A, B, C, D) is a random process vhich is supposed to be 2nd order stationnary in thesense neither its mean and variance depend on A, B, C, D.
  More of this its autocorrelation function is supposed to be a Dirac and its variance finite.
  So E(A, B, C, D) is nothing but a "white noise" W which no longer depend on A, B, C, D.
• In your case you write build Nu from the relation F(A, B, C, D) (your simulated data) but fit the realizations with another model G(A, B, C, D) which differ from F(A, B, C, D) by the fact it contains no quadratic terms.
  So the difference between F(A, B, C, D)  and G(A, B, C, D) is a deterministic function H(A², B², C², D²) and not some random process, and the difference between F(A, B, C, D)  and G_fitted(A, B, C, D) does not come from noisy observations but from a model error (H(A², B², C², D²) ) which is interpreted as noisy observations induced by some hypothetical "stationnary white noise".
• You can of course compute the mean and variance of the residuals  G_fitted(A, B, C, D) - Nu.
  But when the duscrepancy between Nu (= observations) and   G_fitted(A, B, C, D) comes from "a stationnary white noise" of mean, for instance 0, and variance V, one has:
  • Expectation(G_fitted(A, B, C, D) - Nu | A=a, B=b, C=c, D=d)) = 0 whatever a, b, c, d
  • Variance(G_fitted(A, B, C, D) - Nu  |  A=a, B=b, C=c, D=d)) = V whatever a, b, c, d
• Unfortunately this does not apply when Nu = G(A, B, C, D) + H(A², B², C², D²) as in your simulation.
  So usual statistical tests simply should not be used here.
  Specific methods have been developed to handle this situation named as "error model", "error bias", (bayesian framework, model competition, model selection and so on).

I understand your worksheet is more a toy problem than the modeling of a real situation.
But I believe you should generate your observations by writting Nu = G_ß(A, B, C, D) + W and next fit the same model 
G(A, B, C, D). This will give you a fitted model G_ß'(A, B, C, D) where ß' is another set of coefficients which differ from the "true" one ß because of the "noise" W.
About the error you got while running my Maple 2015.2 code: I can'd do any test because I don't have Maple 17.
Nevertheless, the error is induced by Maple 17 not understanding the line

StructuralModel := add(cat(`alpha__`, b) * b, b in basis);

Try replacing it by

StructuralModel := add(alpha[b] * b, b in basis);

Here is the modified worksheet with some add-ons (black text) Nusselt_RSM_Model_mmcdara_2.mw

Let me know if this fixes the error in Maple 17.

@salim-barzani 

I believe there is now a bunch of very skillful Maple users who help you and there is no need for me to intervene, unless on a very specific point or, kike ire, when the answer can be given quickly and that saturdays and sundays calm days on Mapleprimes.

@acer 

Thank you for you return.

@acer 

With Maple 2015.2 solve returns a fonctionion of the form 0 < F And 0 <= F where F is some expression which depends on four parameters (Q_solve.mw): Why does not maple simplify it as 0 < F? (even simplify keeps the result unchanged)?

I noticed this strange result does not appear inyour code.

@justauser 

You're welcome

@Andiguys 

Q_N1_mmcdara.mw

@Andiguys 

Question:
Just a Doubt: Are we optimizing the TR and TP functions? I couldn't find any maximize or NLPsolve  functions in the code you provided.

Answer:
Am I optimizing... ? Of course I do
When you want to maximize some function f(x) it's not mandatoy to use maximize or NLPsolve.
Indeed you can find the set S of solutions of diff(f(x), x)=0 and keep those where eval(diff(f(x), x$2),  x=s) < 0 for all s in S.
If you have an optimization problem with inequality constraints (TR case only), you can do the same thing as before after having replaced f(x) by the lagrangian L(x, a₁, ..., a_K) = f(x) + a₁*g₁(x) + ... + a_K*g_K(x) where K is the number of inequality constraints, g_k(x)  represents the kth inequality constraint and a_k is the kth Lagrange multiplier (you have to compute the gradient of L wrt its K+1 indeterminates and solve a (K+1) by (K+1) algebraic system))

The only thing I didn't do in the TR case is to check which one of the two lolutions retained correspond to a maximum (in my case, and this might depend on the Maple version you use, it's the first solution pairs[1]).
In case of doubts look to the third graph (pairs[1] is blue and pairs[2] is red) to keep the correct solution.

The TP case is that simple that it does not require maximize nor NLPsolve.


By the way: I didn't go into more detail about the Lagragian method in my worksheet because I assumed you were familiar with it from the question you posted a few time ago.

______________________________________________________________________________________


About the discrepancies you point in your question Graph and optimization function giving different...
you're,right, I mede a mistake (plot of diffTRC, Pc) versus Pc instead of TRC vs PC).
Here is the corrected file
Question_New_mmcdara(analytic)_corrected.mw

Here are the values of Pc and TR (in that order) at the maximum for varepsilon=0.09:

[Pc, Max('TRC')] =~ eval(pairs[1], varepsilon=0.09);

         [                                           7]
         [Pc = 1436.173706, Max(TRC) = 6.027683516 10 ]

Sorry for this inconvenience.

______________________________________________________________________________________


About Case 2:
I don't have time to address it right now.
Can you explain simply what are the differences with case ?
If they are marginals, you should be capable to apply yourself the lagrangian method I proposed (which is more suitable than maximize or NLPsolve, IMO).

Then plotting solutions (versus varepsilon) on the same graph is pretty obvious: the only problem is that the different quantities you want to plot have very different orders of magnitude and some of tem will appear as horizontal straight lines.

One workaround is to rescale them, for instance you plot quantity Q1 and, if Q2 is 3 orders of magnitude lower than Q1, you superimpose Q2*1000 instead of Q2 (with legend="Q2*1000")

@one man 

What does "the last ones that you can" mean?
There is an infinity of solutions, os what is for you the "last one"?

For instance

Digits := 25:
fsolve(f, x=10^6+2.*Pi);
eval(f, x=%)
                                              6
                 1.000006283185307179741463 10 
                                          -8
                    2.29978901841373799 10  

You could try and play with the TimeSeriesAnalysis package to detect trends (seasonality), smooth raw data, make predictions and so on.

The "noise" which remains once trend has been removed can be analysed using the SignalProcessing package.

At last the Statistics package could also provide useful functions to refine the analysis.

Interesting work.