% Exchange Algorithm to find the exact optimal design (no replications) % using A-optimal; for mixed 2- and 3-level control and noise factors % All of the 3-level control factors are qualitative variables % Some of the 3-level factors have internal noise (quantitative) %INPUT %kC: kC(1) number of 2-level (qualitative) control factors; kC(2) number of 3-level % (qualitative) control factors %kn: kn number of 2-level noise factors %ki: number of 3-level factor with internal noise %n: number of runs %rho is the correlation parameters defined in Joseph and Delaney (2007) %OUTPUT %D: index of the points selected in the design %Design: returned design %Design=(2-level noise|2-level control|3-level control|3-level factors with %internal noise) %Aopt: the objective value function [D, Design, Aopt]=EAmixedinter(kC,kn,ki,n,rho,wt) p2=kC(1)+kn; % number of 2-level factors p3=kC(2)+ki; % number of 3-level factors N=2^p2*3^p3; %U is the full model matrix %R is the full correlation matrix %points is all the candidates U2=[1,-1;1,1]; % the model matrix for qualitative 3-level factor U3=[1,-(1.5)^0.5,0.5^0.5; 1, 0, -2^0.5; 1, 1.5^0.5, 0.5^0.5]; % the model matrix for quantitative 3-level factor R2=[1,0; 0,(1-rho)/(1+rho)]; % correlation matrix for 3-level qualitative factor R31=[1, 0, 0; 0, (1-rho)/(1+2*rho), 0; 0, 0, (1-rho)/(1+2*rho)]; % correlation matrix for 3-level quantitative factor tempvalue=3+4*rho+2*rho^4; R32=[1, 0, -2^0.5*(rho-rho^4)/tempvalue; 0, 3*(1-rho^4)/tempvalue, 0; -2^0.5*(rho-rho^4)/tempvalue, 0, (rho^4-4*rho+3)/tempvalue]; D2=[-1,1]'; D3=[-1,0,1]'; tempU=1; tempR=1; points=[-1,1]'; for i=1:p2 tempU=kron(U2,tempU); tempR=kron(R2,tempR); end for i=2:p2 points=[kron([1,1]',points),kron(D2,ones(size(points,1),1))]; end for i=1:kC(2) tempU=kron(U3,tempU); tempR=kron(R31,tempR); points=[kron([1,1,1]',points),kron(D3,ones(size(points,1),1))]; end for i=1:ki tempU=kron(U3,tempU); tempR=kron(R32,tempR); points=[kron([1,1,1]',points),kron(D3,ones(size(points,1),1))]; end U=tempU; R=tempR; % log2(flag) is the number of 2-level noise factors contained in the term % log3(flag) is the number of linear-terms generated from any factor with % internal noise contained in the term % log5(flag) is the number of quadratic-terms generated from any factor % with internal noise contained in the term flag=1; for i=1:kn flag=kron([1,2],flag); end for i=1:kC(1) flag=kron([1,1],flag); end for i=1:kC(2) flag=kron([1,1,1],flag); end for i=1:ki flag=kron([1,3,5],flag); end %pick up the initial points. n0=floor(n/2); pick=0; D=myrandint(1,n0,1:N,'noreplace'); % introduce some randomness to avoid trapping in the local optimal temp=U(D,:)*R; V=temp*U(D,:)'; % as long as D has no replicates, V is not singular M=temp'*inv(V)*temp; opt=wt*diag(M); left=N-n0; candidate=setdiff(1:N,D); % all the left candidate points %choose the next n-n0 points l=n0+1; while l<=n F=R-M; delta0=0; i=1; while i<=left f=U(candidate(i),:); a=f*F; d=a*f'; delta1=wt*(a.*a)'/d; if delta1>delta0 a0=a; d0=d; pick=candidate(i); delta0=delta1; end i=i+1; end M=M+(a0'*a0)/d0; D=[D,pick]; candidate=candidate(candidate~=pick); l=l+1; opt=opt+delta0; left=left-1; end while 1 compopt=opt; for i=1:n pick=0; D1=setdiff(D,D(i)); temp=U(D1,:)*R; V=inv(temp*U(D1,:)'); M1=temp'*V*temp; opt1=wt*diag(M1); F=R-M1; delta0=opt-opt1; for j=1:left f=U(candidate(j),:); a=f*F; d=a*f'; delta1=wt*(a.*a)'/d; if delta1>delta0 pick=candidate(j); delta0=delta1; end end if pick~=0 opt=opt1+delta0; candidate=candidate(candidate~=pick); candidate=[candidate,D(i)]; D(i)=pick; end end if abs(compopt-opt)<1e-8 break end end D=sort(D); Aopt=opt/(wt*diag(R)); Design=points(D,:);