MATLAB: How to determine the stationary points of a fitted second order model using Matlab

code generationodeoptimization

Below is the Matlab code developed for fitting a second order model to the data X1,X2,X3 and X4
%Matlab algorithm for Fitting the Second Order Model for the Stress at the Manifold X1=[3688,5229,3688,5229,3688,5229,3688,5229,3688,5229,3688,5229,3688,5229,3688,5229 ,2918,6000,4459,4459,4459,4459,4459,4459,4459,4459]'; X2=[2761,2761,3587,3587,2761,2761,3587,3587,2761,2761,3587,3587,2761,2761,3587,3587,3174,3174,2348,4000,3174,3174,3174,3174,3174,3174]'; X3=[3761,3761,3761,3761,6587,6587,6587,6587,3761,3761,3761,3761,6587,6587,6587,6587, 5174,5174,5174,5174,2348,8000,5174,5174,5174,5174]'; X4=[2780,2780,2780,2780,2780,2780,2780,2780,3646,3646,3646,3646,3646,3646,3646,3646, 3213,3213,3213,3213,3213,3213,2348,4078,3213,3213]'; Stress=[246,255,248,257,243,243,245,246,245,257,249,261,243,244,245,248, 247,270,244,248,249,243,246,247,247,247]'; % This algorithm creates the design Matrix X=[X1 X2 X3 X4]; % Algorithm for creating the second order quadratic Model for the Stress at the Manifold format short e; mdl=LinearModel.fit(X, Stress, 'quadratic');
% Response surface plots(3-D) showing the effect of the variation of the run lenghts of the Jumper on the Stress at the Manifold for the first order model %This algorithm plots X1 and X2 against the response while X3 and X4 and fixed at zero [x1,x2] = meshgrid(linspace(2918,6000,15), linspace(2348,4000,15)); A = zeros(225,4); A(:,1) = x1(:); A(:,2) = x2(:); B = predict(mdl,A); B = reshape(B,size(x1)); %This algorithm plots X1 and X3 against the response while X2 and X4 and fixed at zero [x1,x3] = meshgrid(linspace(2918,6000,15), linspace(2348,8000,15)); C = zeros(225,4); C(:,1) = x1(:); C(:,2) = x3(:); D = predict(mdl,C); D = reshape(D,size(x1));
%This algorithm plots X1 and X4 against the response while X2 and X3 and fixed at zero [x1,x4] = meshgrid(linspace(2918,6000,15), linspace(2348,4078,15)); E = zeros(225,4); E(:,1) = x1(:); E(:,2) = x4(:); F = predict(mdl,E); F = reshape(F,size(x1));
%This algorithm plots X2 and X3 against the response while X1 and X4 and fixed at zero [x2,x3] = meshgrid(linspace(2348,4000,15), linspace(2348,8000,15)); G = zeros(225,4); G(:,1) = x2(:); G(:,2) = x3(:); H = predict(mdl,G); H = reshape(H,size(x2));
%This algorithm plots X2 and X4 against the response while X1 and X3 and fixed at zero [x2,x4] = meshgrid(linspace(2348,4000,15), linspace(2348,4078,15)); I = zeros(225,4); I(:,1) = x2(:); I(:,2) = x4(:); J = predict(mdl,I); J = reshape(J,size(x2));
%This algorithm plots X3 and X4 against the response while X1 and X2 and fixed at zero [x3,x4] = meshgrid(linspace(2348,8000,15), linspace(2348,4078,15)); K = zeros(225,4); K(:,1) = x3(:); K(:,2) = x4(:); L = predict(mdl,K); L = reshape(L,size(x3));
subplot(2,3,1); surf(x1,x2,B),view(-150,20), xlabel('X1'), ylabel ('X2'), zlabel ('Stress'),title('(a) Fixed levels: x3=0 x4=0'); subplot(2,3,2); surf(x1,x3,D), view(-150,20), xlabel('X1'), ylabel ('X3'), zlabel ('Stress') ,title('(b) Fixed levels: x2=0 x4=0'); subplot(2,3,3); surf(x1,x4,F), view(-150,20), xlabel('X1'), ylabel ('X4'), zlabel ('Stress') ,title('(c) Fixed levels: x2=0 x3=0'); subplot(2,3,4); surf(x2,x3,H), view(-150,20), xlabel('X2'), ylabel ('X3'), zlabel ('Stress') ,title('(d) Fixed levels: x1=0 x4=0'); subplot(2,3,5); surf(x2,x4,J), view(-150,20), xlabel('X2'), ylabel ('X4'), zlabel ('Stress') ,title('(e) Fixed levels: x1=0 x3=0'); subplot(2,3,6); surf(x3,x4,L), view(-150,20), xlabel('X3'), ylabel ('X4'), zlabel ('Stress') ,title('(f) Fixed levels: x2=0 x3=0');
It is required to determine the optimum settings of the variables(X1,X2,X3 and X4) that will the minimum response(Stress)

Best Answer

Without a bound on the variables X1 through X4, there is no minimum of the model response. The fitted model has negative eigenvalues.
You can see this by examining the coefficients:
 b = mdl.Coefficients.Estimate
Place the quadratic terms in a matrix and examine the eigenvalues.
 b3 = diag(b(12:15));
 b3(1,2:4) = b(6:8);
 b3(2,3:4) = b(9:10);
 b3(3,4) = b(11);
 b3 = b3 + b3';
 eig(b3)
 ans =
   1.0e-05 *
   -0.5606
   -0.1965
   -0.0831
    0.9961
If I made a mistake, and counted the diagonal terms twice, well, redo the calculation with half the diagonal, and the answer still has negative eigenvalues.
Alan Weiss
MATLAB mathematical toolbox documentation
Related Question