Subscribe Bookmark RSS Feed

ODE for high order response surface model

stochastic_mind

Community Trekker

Joined:

Jul 1, 2015

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

SM

13 REPLIES
ian_jmp

Staff

Joined:

Jun 23, 2011

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

stochastic_mind

Community Trekker

Joined:

Jul 1, 2015

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,

SM

Peter_Bartell

Joined:

Jun 5, 2014

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

Community Trekker

Joined:

Jul 1, 2015

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,

SM

billw_jmp

Staff

Joined:

Jul 2, 2014

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

stochastic_mind

Community Trekker

Joined:

Jul 1, 2015

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,

SM

Peter_Bartell

Joined:

Jun 5, 2014

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

Community Trekker

Joined:

Jul 1, 2015

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,

SM

Peter_Bartell

Joined:

Jun 5, 2014

I just asked JMP's Custom Design platform to find me a design with 10 continuous factors, a model of all main effects,  all 2 way interactions, and 3 way interactions and it took about 15 seconds to converge to a design. So I think you should be OK? Do you have JMP Pro? With all these disparate models you'll be fitting, you might want to take a look at all of the results within the framework of the Model Comparison platform which is available in JMP Pro only. You can do lots of design diagnostics within the DOE platform too...but if you want to take it all the way through modeling with actual results, the Model Comparison platform might be very efficient for you?