nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

ODE for high order response surface model

stochastic_mind

Dear all:

I need to construct classical optimal designs (A-, D-, E- and G-optimality criteria) for a high order response surface model with continuous factors. My response surface model is a multidimensional polynomial and it has about 300 regression parameters to be determined.

My understanding is that exchange methods for this problem are not useful since the search process is very time consuming.

I would appreciate if you could let me know about the best algorithm for this problem to construct optimal designs efficiently. Is JMP capable of constructing optimal designs for my problem?

Best

Accepted answers

All comments

ian_jmp

If you know the model you want to fit, JMP Custom Designer will help with the D-optimal (or I-optimal) case: Custom Designs

Peter_Bartell

I guess I'm more curious than anything else...is your number of parameter estimates driven by lots and lots of factors with relatively low order interaction/polynomial terms or just a smaller number of factors but a desire to estimate very high order interaction effects? Sounds like you have a system where the usual conditions around effect sparsity don't apply? If nothing else sounds like an interesting system under study that would necessitate either sets of experimental requirements.

stochastic_mind

Ian:

Thanks for your reply. It seems that the highest power is 5 which would not be large enough for my model. Other than that, is there any way that I can enter my response model manually similar to what they do in Gosset (http://neilsloane.com/gosset/)?

Best,

stochastic_mind

Peter:

Great question! Please let me be more specific about my system. In fact, I have the both scenarios you mentioned. In one case, I have a large number of continuous factors with low order polynomials. Say I have 10 factors; x_1, …, x_10 \in [-1,1]. Let's say least square regression problem is given by:

u(x_1, …, x_10) = \sum_{i=1}^{i=P} a_i*y_i(x_1, …, x_10)

where u is the response, y_i is the basis function and a_i's are the regression coefficients to be determined. P is also the number of unknown regression coefficients. The second order y_i is given by:

y_i = 0.0541266*x_1 + 0.0541266*x_2 + 0.0541266*x_3 + 0.0541266*x_4 + 0.0541266*x_5 + 0.0541266*x_6 + 0.0541266*x_7 + 0.0541266*x_8 + 0.0541266*x_9 + 0.0541266*x_10 + 0.09375*x_1*x_2 + 0.09375*x_1*x_3 + 0.09375*x_1*x_4 + 0.09375*x_2*x_3 + 0.09375*x_1*x_5 + 0.09375*x_2*x_4 + 0.09375*x_1*x_6 + 0.09375*x_2*x_5 + 0.09375*x_3*x_4 + 0.09375*x_1*x_7 + 0.09375*x_2*x_6 + 0.09375*x_3*x_5 + 0.09375*x_1*x_8 + 0.09375*x_2*x_7 + 0.09375*x_3*x_6 + 0.09375*x_4*x_5 + 0.09375*x_1*x_9 + 0.09375*x_2*x_8 + 0.09375*x_3*x_7 + 0.09375*x_4*x_6 + 0.09375*x_1*x_10 + 0.09375*x_2*x_9 + 0.09375*x_3*x_8 + 0.09375*x_4*x_7 + 0.09375*x_5*x_6 + 0.09375*x_2*x_10 + 0.09375*x_3*x_9 + 0.09375*x_4*x_8 + 0.09375*x_5*x_7 + 0.09375*x_3*x_10 + 0.09375*x_4*x_9 + 0.09375*x_5*x_8 + 0.09375*x_6*x_7 + 0.09375*x_4*x_10 + 0.09375*x_5*x_9 + 0.09375*x_6*x_8 + 0.09375*x_5*x_10 + 0.09375*x_6*x_9 + 0.09375*x_7*x_8 + 0.09375*x_6*x_10 + 0.09375*x_7*x_9 + 0.09375*x_7*x_10 + 0.09375*x_8*x_9 + 0.09375*x_8*x_10 + 0.09375*x_9*x_10 + 0.104816*x_1^2 + 0.104816*x_2^2 + 0.104816*x_3^2 + 0.104816*x_4^2 + 0.104816*x_5^2 + 0.104816*x_6^2 + 0.104816*x_7^2 + 0.104816*x_8^2 + 0.104816*x_9^2 + 0.104816*x_10^2 - 0.318136

For the second case, I have only two factors x_1 and x_2 \in [-1,1] with a high order polynomial, i.e., y_i is given by:

y_i = 3.11846*x_1 + 3.11846*x_2 + 5.5*(1.875*x_1 - 8.75*x_1^3 + 7.875*x_1^5)*(1.875*x_2 - 8.75*x_2^3 + 7.875*x_2^5) + 3.3541*(1.5*x_2^2 - 0.5)*(4.375*x_1^4 - 3.75*x_1^2 + 0.375) + 3.3541*(1.5*x_1^2 - 0.5)*(4.375*x_2^4 - 3.75*x_2^2 + 0.375) + 3.1225*x_2*(6.5625*x_1^2 - 19.6875*x_1^4 + 14.4375*x_1^6 - 0.3125) + 3.1225*x_1*(6.5625*x_2^2 - 19.6875*x_2^4 + 14.4375*x_2^6 - 0.3125) - 4.7697*(1.5*x_2 - 2.5*x_2^3)*(6.5625*x_1^2 - 19.6875*x_1^4 + 14.4375*x_1^6 - 0.3125) - 4.7697*(1.5*x_1 - 2.5*x_1^3)*(6.5625*x_2^2 - 19.6875*x_2^4 + 14.4375*x_2^6 - 0.3125) + 2.87228*x_2*(1.875*x_1 - 8.75*x_1^3 + 7.875*x_1^5) + 2.87228*x_1*(1.875*x_2 - 8.75*x_2^3 + 7.875*x_2^5) + 3.77492*x_2*(2.46094*x_1 - 36.0938*x_1^3 + 140.766*x_1^5 - 201.094*x_1^7 + 94.9609*x_1^9) + 3.77492*x_1*(2.46094*x_2 - 36.0938*x_2^3 + 140.766*x_2^5 - 201.094*x_2^7 + 94.9609*x_2^9) - 4.38748*(1.5*x_2 - 2.5*x_2^3)*(1.875*x_1 - 8.75*x_1^3 + 7.875*x_1^5) - 4.38748*(1.5*x_1 - 2.5*x_1^3)*(1.875*x_2 - 8.75*x_2^3 + 7.875*x_2^5) + 1.5*x_1*x_2 - 2.29129*x_2*(1.5*x_1 - 2.5*x_1^3) - 2.29129*x_1*(1.5*x_2 - 2.5*x_2^3) + 4.5*(4.375*x_1^4 - 3.75*x_1^2 + 0.375)*(4.375*x_2^4 - 3.75*x_2^2 + 0.375) + 4.03113*(1.5*x_2^2 - 0.5)*(6.5625*x_1^2 - 19.6875*x_1^4 + 14.4375*x_1^6 - 0.3125) + 4.03113*(1.5*x_1^2 - 0.5)*(6.5625*x_2^2 - 19.6875*x_2^4 + 14.4375*x_2^6 - 0.3125) + 3.5*(1.5*x_1 - 2.5*x_1^3)*(1.5*x_2 - 2.5*x_2^3) + 3.7081*(1.5*x_2^2 - 0.5)*(1.875*x_1 - 8.75*x_1^3 + 7.875*x_1^5) + 3.7081*(1.5*x_1^2 - 0.5)*(1.875*x_2 - 8.75*x_2^3 + 7.875*x_2^5) + 3.57071*x_2*(54.1406*x_1^4 - 9.84375*x_1^2 - 93.8438*x_1^6 + 50.2734*x_1^8 + 0.273438) + 3.57071*x_1*(54.1406*x_2^4 - 9.84375*x_2^2 - 93.8438*x_2^6 + 50.2734*x_2^8 + 0.273438) + 1.93649*x_1*(1.5*x_2^2 - 0.5) + 1.93649*x_2*(1.5*x_1^2 - 0.5) - 2.95804*(1.5*x_1 - 2.5*x_1^3)*(1.5*x_2^2 - 0.5) - 2.95804*(1.5*x_2 - 2.5*x_2^3)*(1.5*x_1^2 - 0.5) - 3.3541*x_2*(2.1875*x_1 - 19.6875*x_1^3 + 43.3125*x_1^5 - 26.8125*x_1^7) - 3.3541*x_1*(2.1875*x_2 - 19.6875*x_2^3 + 43.3125*x_2^5 - 26.8125*x_2^7) + 5.12348*(1.5*x_2 - 2.5*x_2^3)*(2.1875*x_1 - 19.6875*x_1^3 + 43.3125*x_1^5 - 26.8125*x_1^7) + 5.12348*(1.5*x_1 - 2.5*x_1^3)*(2.1875*x_2 - 19.6875*x_2^3 + 43.3125*x_2^5 - 26.8125*x_2^7) + 5.40833*(4.375*x_2^4 - 3.75*x_2^2 + 0.375)*(6.5625*x_1^2 - 19.6875*x_1^4 + 14.4375*x_1^6 - 0.3125) + 5.40833*(4.375*x_1^4 - 3.75*x_1^2 + 0.375)*(6.5625*x_2^2 - 19.6875*x_2^4 + 14.4375*x_2^6 - 0.3125) + 4.97494*(1.875*x_1 - 8.75*x_1^3 + 7.875*x_1^5)*(4.375*x_2^4 - 3.75*x_2^2 + 0.375) + 4.97494*(1.875*x_2 - 8.75*x_2^3 + 7.875*x_2^5)*(4.375*x_1^4 - 3.75*x_1^2 + 0.375) + 4.60977*(1.5*x_2^2 - 0.5)*(54.1406*x_1^4 - 9.84375*x_1^2 - 93.8438*x_1^6 + 50.2734*x_1^8 + 0.273438) + 4.60977*(1.5*x_1^2 - 0.5)*(54.1406*x_2^4 - 9.84375*x_2^2 - 93.8438*x_2^6 + 50.2734*x_2^8 + 0.273438) + 2.5*(1.5*x_1^2 - 0.5)*(1.5*x_2^2 - 0.5) + 2.59808*x_2*(4.375*x_1^4 - 3.75*x_1^2 + 0.375) + 2.59808*x_1*(4.375*x_2^4 - 3.75*x_2^2 + 0.375) - 4.33013*(1.5*x_2^2 - 0.5)*(2.1875*x_1 - 19.6875*x_1^3 + 43.3125*x_1^5 - 26.8125*x_1^7) - 4.33013*(1.5*x_1^2 - 0.5)*(2.1875*x_2 - 19.6875*x_2^3 + 43.3125*x_2^5 - 26.8125*x_2^7) + 18.6023*x_1^2 - 51.7429*x_1^3 + 18.6023*x_2^2 - 186.095*x_1^4 - 51.7429*x_2^3 + 235.976*x_1^5 - 186.095*x_2^4 + 638.9*x_1^6 + 235.976*x_2^5 - 386.351*x_1^7 + 638.9*x_2^6 - 875.481*x_1^8 - 386.351*x_2^7 + 206.963*x_1^9 - 875.481*x_2^8 + 413.407*x_1^10 + 206.963*x_2^9 + 413.407*x_2^10 - 3.96863*(1.5*x_1 - 2.5*x_1^3)*(4.375*x_2^4 - 3.75*x_2^2 + 0.375) - 3.96863*(1.5*x_2 - 2.5*x_2^3)*(4.375*x_1^4 - 3.75*x_1^2 + 0.375) - 0.6201

I know that this is not a standard approach in optimal design of experiments. In fact, my system is described by a computer simulation rather than a real world experiment. That is why I need to have higher order terms.

Is it even computationally feasible to construct an optimal design for such high order polynomials?

Best,

Bill_Worley

SM,

Since you are working with computer simulations have you considered Space Filling designs? With a response surface design you will likely get multiple center point conditions and with computer simulations you only need one point per condition due to the deterministic nature of the simulation. That is unless you are forcing variation. From past experience computer simulations are highly non-linear and using a space filing design will help you better understand and ultimately model that non-linearity.

Gaussian Process and Neural Net models are typically better suited to handle the modeling for computer simulations to build your surrogate model. Neural Nets are typically faster than Gaussian Process especially if you are up past 100 simulation conditions. RS models can be nice if they can fit your data because they are usually simpler than either Gaussian or NN.

The Fast Flexible Filling design with MaxPro optimality will give you a good option versus a response surface design especially for a design with higher dimensions.

Best,

Bill

Peter_Bartell

Ahhh...the computer simulation scenario makes sense to me. I was wondering even with a relatively small number of factors, how in the world are you going to keep track of factor settings, run order, keeping responses tied to actual experimental outcome units, and the 1001 logistical details that can derail even, say a 2**k-p, 32 run design for the required degrees of freedom for a 300+ parameter estimate model. Just too big a logistical challenge to even consider.

stochastic_mind

Peter:

The model with 300 regression coefficients would be my extreme case. The model I provided (y_i) has 66 unknowns. What about this model? Is it still too large and challenging? I have to add that I need to find optimal deigns with a number of runs from one to ten times larger than the number of unknowns.

Best,

stochastic_mind

Bill:

Thank you for your reply. I have tested Latin Hypercube designs and they work well for my problem. In fact, I am not just trying to find an optimal design for my problem, I want to see how different optimality criteria affect the accuracy of the solution and how they are compared for example with LH designs.

I am a newbie in ODE, so, my major question would be: is it even feasible to construct D-, A-, E, and G-optimal designs for the model described in reply to Peter's answer?

Best,

Mark_Bailey

The various space filling design options assume that the computer simulation is entirely deterministic, so most criteria that might be used for an optimal design to collect data from a stochastic system do not apply. For example, D, A, G, and I do not apply because they aim to reduce the standard error of your estimates of the model parameters or mean response. You have standard errors of zero.

Your game is interpolation, so space filling designs are the way to go and the Gaussian process and neural network models are the best interpolators. Forget about using very high order polynomials and linear regression.

stochastic_mind

Mark:

Thank you for your email. My computer simulation is a stochastic one, so the standard error would not be zero. The uniform factors I have in my response model are basically random variables.

I am trying to figure out if there is an algorithm/software which can construct alphabetic optimal designs for high order polynomials. Are you aware of such an algorithm? I know that exchange methods are not suitable for my model since it is a large scale problem and those methods are greedy and time consuming. But what about simulated annealing or genetic algorithms? Are they capable of constructing optimal designs for high order polynomials (say up to 10^th order)?

Best,

bradleyjones

SM - It is possible to construct designs for higher order polynomials using JMP's Custom Design tool and scripting. Here is a sample script.

DOE(Custom Design,{Add Response( Maximize, "Y", ., ., . ),

Add Factor( Continuous, -1, 1, "X1", 0 ), Set Random Seed( 1318892398 ),

Add Term( {1, 0} ), Add Term( {1, 1} ), Add Term( {1, 2} ), Add Term( {1, 3} ),

Add Term( {1, 4} ), Add Term( {1, 5} ), Add Term( {1, 6}),Set Sample Size( 24 )}

);

Note that the last Add Term command adds a sixth order polynomial term. The current maximum order for any term is 10 but the UI only surfaces terms up to 5th order.