% Exchange Algorithm to find the exact optimal design (no replications) % For 2-level control and 2-level noise factors only % INPUT % kC: number of control factors; % kn: number of noise variables % n: number of runs % r: parameter of the prior correlation matrix of beta (usually set to be % 1/3) % OUTPUT % D: index of the points selected in the design % Design: design matrix=(control factors | noise factors) % Aopt: the objective value function [D,Design, Aopt]=EA2level(kC,kn,n,r) p=kC+kn; %p is the number of factors N=2^p; %N is the number of candidates %U is the full model matrix for all the candidate points %points is the design matrix (all the candidates points) %R is the correlation matrix %points is the design settings for all the factors. U=[1,-1;1,1]; R=diag([1,r]); points=[-1,1]'; for i=2:p U=[U,-U;U,U]; points=[kron([1,1]',points),kron([-1,1]',ones(size(points,1),1))]; temp=diag(R)'; R=diag([temp,r*temp]); end %flag is used to indicate how many noise factors in the term =log2(flag) wt=zeros(1,N); flag=1; for i=1:kC flag=kron([1,1],flag); end for i=1:kn flag=kron([1,2],flag); end % wt is the diagonal matrix A, wt is user specified wt(flag==2)=1; % pick up the initial points n0=floor(n/3); 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,:);