% %SEP2.m % % function [choice,maxnac]=sep2(Z,isdiscrete,x0,f0,maxnac0) % % optimization via seperable underestimation % for min f(x), s.t. x\in T % where T=S1xS2, S1 a finite subset of R^{n_d}, S2 a box in R^{n_c}. % % INPUTS: % % Z is a cellarray of tables, represents T, constrains the possible choices, i.e., % Z{i}=[a b], an interval in the continuous case % Z{i} has size N_i x n_i, n_i the # of design variables, N_i the # of entries in the table Z{i}, % in the discrete case % % isdiscrete is a binary 1 x length(Z) vector with entry 1 iff variable i is discrete % % x0 are evaluation points, should have size Nxn, % where N is the # of evaluations, n the length of Z, i.e., n=n_d+n_c % % f0 are the function values at x0, should have size 1xN % % optional input: maxnac0, the max number of active constraints in the algorithm, cf. [1] % % OUTPUTS: % % choice is the suggested optimal solution % % optional output: maxnac, the max number of active constraints in the algorithm % % REQUIREMENTS: % % CVX: Matlab Software for Disciplined Convex Programming % http://www.stanford.edu/~boyd/cvx/ % % BIBLIOGRAPHY: % % [1] M. Fuchs and A. Neumaier. Autonomous robust design optimization with potential clouds. % International Journal of Reliability and Safety, % Special Issue on Reliable Engineering Computing, 2008. % [2] M. Fuchs. Uncertainty modeling in higher dimensions: Towards robust design optimization. % PhD thesis, University of Vienna, Faculty of Mathematics, 2008. % % DOWNLOAD SOURCE: % % martin-fuchs.net function [choice,maxnac]=sep2(Z,isdiscrete,x0,f0,maxnac0) if nargin==3 maxnac0=inf; end; choice=[]; maxnac=[]; % unique x0str=[]; for i=1:size(x0,1) x0str{i}=num2str(x0(i,:)); end; [temp,per]=unique(x0str); % sort, unique x0=x0(per,:); f0=f0(per); lf0=length(f0); if lf0==0 || isempty(x0) return; end; lZ=length(Z); for i=1:lZ lZZ(i)=size(Z{i},1); end; % find index of Z with interval design vars, and max table length contind=[]; maxind=0; for i=1:lZ if ~isdiscrete(i) contind=[contind i]; end; if size(Z{i},1)>maxind maxind=size(Z{i},1); end; end; intcheck=(length(contind)==lZ); % is 1 if all choices are continous infi=max(f0); % inf for infeasible choices in cvx f=f0; for i=1:lf0 % scale f0 if isnan(f0(i)) || isinf(f0(i)) f(i)=1; else f(i)=f0(i)/infi; end; end; % transform the choices x0 to constraint lhs constr=cell(0,0); for i=1:lf0 constr{i}=''; for j=1:lZ if any(j==contind) constr{i}=[constr{i} num2str(x0(i,j)) '*(q' num2str(j) '(1)+' num2str(x0(i,j)) '*q' num2str(j) '(2))+']; else constr{i}=[constr{i} 'q' num2str(j) '(' num2str(x0(i,j)) ')+']; end; end; if intcheck constr{i}=[constr{i} 'c+']; end; constr{i}=constr{i}(1:end-1); end; % objective function string objstr=[]; for i=1:lZ objstr=[objstr 'q' num2str(i) '''*q' num2str(i) '+'];%[objstr 'sum(qa' num2str(i) ')+']; end; if intcheck objstr=[objstr 'c^2+']; end; objstr=objstr(1:end-1); % max number of active constraints to avoid overestimation/infeasibility/overdetermination % in the cvx problem: if maxnac0==0 maxnac=inf; else maxnac=min(min(lZ+sum(lZZ)-1,maxnac0),lf0); end; % index of constraints permuted ind=1:lf0; lindtemp=min(maxnac,lf0); % take equidistant constr numbers if lindtemp==1 indtemp=ind(1); else indtemp=ind(ceil([1:(lf0-1)/(lindtemp-1):lf0])); end; lind0=0; firstiter=1; niter=1; cvxinf=0; while 1 % loop until all constraints <= f0 if firstiter==1 while 1 % loop to get the maximum<=macnac of active constraints for one underestimator cvx_begin cvx_quiet(true); for i=1:lZ if any(i==contind) variable(['q' num2str(i) '(2)']); else variable(['q' num2str(i) '(' num2str(lZZ(i)) ')']); end; end; if intcheck warning('No integer choice variables found.'); variable('c(1)'); end; minimize(eval(objstr)); subject to for i=1:lindtemp eval(constr{indtemp(i)})==f(indtemp(i)); end; cvx_end %cvx_status if isequal(cvx_status,'Infeasible') || isequal(cvx_status,'Failed') || isequal(cvx_status,'Overdetermined') %|| isequal(cvx_status,'Inaccurate/Solved') %|| isempty(ind1) lindtemp=floor((lind0+lindtemp)/2); % take equidistant constr numbers if lindtemp==1 indtemp=ind(1); else indtemp=ind(ceil([1:(lf0-1)/(lindtemp-1):lf0])); end; else lind0=length(indtemp); lindtemp=floor((lind0+lindtemp)/2); % take equidistant constr numbers if lindtemp==1 indtemp=ind(1); else indtemp=ind(ceil([1:(lf0-1)/(lindtemp-1):lf0])); end; for i=1:lZ eval(['qq' num2str(i) '=q' num2str(i) ';']); % last feasible solution end; end; % stopping criterion if lindtemp==lind0 maxnac=lind0; break; end; end; lindtemp=maxnac; if maxnac==0 choice=[]; return; end; for i=1:lZ eval(['q' num2str(i) '=qq' num2str(i) ';']); % from qq back to q end; else cvx_begin cvx_quiet(true); for i=1:lZ if any(i==contind) variable(['q' num2str(i) '(2)']); else variable(['q' num2str(i) '(' num2str(lZZ(i)) ')']); end; end; if intcheck warning('No integer choice variables found.'); variable('c(1)'); end; minimize(eval(objstr)+eval(viopen)); subject to for i=1:lindtemp eval(constr{indtemp(i)})==f(indtemp(i)); end; cvx_end %cvx_status if isequal(cvx_status,'Infeasible') || isequal(cvx_status,'Failed') || isequal(cvx_status,'Overdetermined') %|| isequal(cvx_status,'Inaccurate/Solved') %|| isempty(ind1) for i=1:lZ eval(['q' num2str(i) '=qq' num2str(i) ';']); % from qq back to q end; cvxinf=1; else for i=1:lZ eval(['qq' num2str(i) '=q' num2str(i) ';']); % from qq back to q end; cvxinf=0; end; end; if cvxinf==0 % compute which constraints have been violated the most (maxnac many) vio=zeros(1,lf0); for i=ind vio(i)=f(i)-eval(constr{i}); end; [temp,per]=sort(vio); indtemp=per(1:maxnac); % violation penalty function for the objective function %% indviopen=per(lf0-maxnac+1:lf0); viopen=[]; for i=indviopen viopen=[viopen '(' constr{i} '-f(' num2str(i) '))^2+']; end; viopen=viopen(1:end-1); else % take a subset of constraints indtemp=per(1:ceil(lindtemp/2)); lindtemp=length(indtemp); end; % stopping criterion (all underestimating or too many iterations) if temp(1)>=-eps if firstiter==0 || temp(end)<=1e-6 % cvx precision break; end; elseif niter>=lf0/maxnac && maxnac