function [aMix] = soptim(aMix, aMixu, ufc, nstep, chis) % soptim performs simultaneous advisory design for normal mixture % % [aMix] = soptim(aMix, aMixu, ufc, nstep, chis) % [aMix] = soptim(aMix, aMixu, ufc, nstep) chis = 1 % [aMix] = soptim(aMix, aMixu, ufc) nstep = [200, 1] % % aMix : advised mixture of the type ARX LS enriched on following control states: % strc : common control structure % ufc : normalised vector qualifying components: % dangerous component (0), not dangerous (positive number) % kc : lift of quadratic forms % UDc : cell vector of u'du decompositions of KLD kernels % udca : u'du decomposition of average KLD kernel in UDc % kca : average lift of quadratic forms kc % aMixu : desired mixture (user's target) of the type ARX LS with control states % ufc : vector qualifying components: 0 - dangerous component, (1) - not % nstep : parameters [ns1,per] determining design horizon, i.e. horizon = ns1*per; % ns1 : number of block repetition % per : horizon of a block % if nstep is defined by parameter nsl only then per is set to 1 % chis : indicates strategy chosen: chis=1 for receding horizon (default) and chis=-1 for IST % % Design : J. Bohm % Updated : June, 2002 % Project : ProDaCTools, IST-1999-12058 % See also : udupdt, getdvect, facchng, facarxls % References : \ref{ch9} % Note : % Updated : if (nargin < 3) | ~any(ufc)error(sprintf('%s\n%s\n%s','ufc was not correcly set, define it as a vector','of a length "ncom" having at least one nonzero element ','or use a function ufcgen')) ;end if (nargin <4), nstep=[200,1];chis=1;end if (nargin <5), chis=1;end if(isempty(aMix.states.uchn)) error('uoptim needs nonempty list of channels with recognisable actions'); end; % normalisation of ufc ufc = ufc/sum(ufc); %Inititialization ncom = length(aMix.dfcs); dfcs = aMix.dfcs; strc = aMix.states.strc; % common control structure nPsi = max(size(strc)); % length of regression vector + data pochn = aMix.states.pochn; % list of channels with o-innovations npochn = length(pochn); % number of channels with o-innovations nychn = length(aMix.states.modelled); % number of modelled channels nouts = length(aMix.states.outs); % number of innovation channels npsi = nPsi-nychn; % length of the regression vector kc0 = aMixu.states.kc; udca = aMixu.states.udca; kca = aMixu.states.kca; coms = aMix.coms; Ethz = zeros(1,nPsi); UDc = aMixu.states.UDc; lss = length(nstep); lrica = zeros(1,nPsi-1); %test of aMixu coves=zeros(1,npochn); % if npochn~=length(aMixu.Facs),error('aMixu not correctly set'); end for i=1:npochn, coves(i)=aMixu.Facs{i}.cove; end if ~any(coves), error('aMixu not correctly set, Facs{.}.cove must be >0'),end for i=1:ncom, % cycle over number of components ncom Ric{i}= zeros(nPsi); % KLD kernels kcc(i)=0; % lift of quadratic forms end df=dfcs/sum(dfcs); % setting of design horizon if lss==2 steps = nstep(1)*nstep(2); per=nstep(2); else steps = nstep; per=1; end if chis>0, % if the strategy starts from zero udca = zeros(nPsi); kca =0; end %Main design cycle, iterations over the horizon of the criterion for iter=1:steps, % ===================== iterations till design horizon % ricmn is an auxiliary array accumulating results of optimization if mod(iter-1,per)==0, % shift of a matrix from bottom right to top left by nychn % if nPsi>nychn+1 [udca, lrica]= ricshift(udca,lrica,nychn,nPsi,npsi); % end %if nPsi udca(nPsi,nPsi)= 0; end %if mod for i=1:ncom % ......................... cycle over all components if mod(iter-1,per)==0, % if ufc(i)==0, kcc(i)=1e30; continue; end % excluding bad components ric = udca; % ric is auxiliarry working array lric(i,:)=lrica; kcc(i) = -npochn +kc0(i); % adding to each component its stationary loss for j=1:npsi, % red=UDc{i}(j,:); red(j)=1; ric= udupdt(ric,red,ufc(i)*UDc{i}(j,j)); end % for j else % iteration continues in corresponding component kernel and lift % shift of a matrix from bottom right to top left by nychn if nPsi>nychn+1, [ric,lric(i,:)]=ricshift(Ric{i},lric(i,:),nychn,nPsi,npsi); end % if nPsi ric(nPsi,nPsi) = Ric{i}(nPsi,nPsi); end %if mod % expectation is calculated channel by channel for j=1: nouts %--------------------- cycle over innovation channels indv = ~isempty(find(strc(1,j)==pochn)); % indicator if the channel is o-innovation [ric,lric,kcc]= ricexp(ric,lric,kcc,i,j,aMix,nPsi); % the penalization is used if indv, % visibility indicator [ric,lric,kcc]= ricpen(ric,lric,kcc,i,j,aMix,aMixu,nPsi); end %end if indv end % ------------------------------ reduced all factors of i-th component % now penalization for j=nouts+1: nychn, [ric,lric,kcc]= ricpenu(ric,lric,kcc,i,j,aMix,aMixu,nPsi); % if ric(j,j)>eps, % kcc(i)=kcc(i)-lric(i,j)*lric(i,j)/ric(j,j)/4; % disp('pred m'); % keyboard % % lric(i,:)=lric(i,:)-lric(i,j)*[zeros(1,j) ric(j,j+1:end-1)]; % disp('po m'); % keyboard % end % if end %end for [l,d]= ld2ld(ric(nouts+1: nychn,nouts+1: nychn)); lric(i,nychn+1:end)=lric(i,nychn+1:end)-lric(i,nouts+1: nychn)*inv(l)*ric(nouts+1: nychn,nychn+1:end-1); Ric{i} = ric; end % ...................... done for all components % now putting losses together if mod(iter,per)==0, udca=zeros(size(ric)); lrica = zeros(1,nPsi-1); det=ones(1,ncom); %%putting it together for i=1:ncom, % if ufc(i)==0, continue; end for j=nouts+1:nPsi, red = Ric{i}(j,:); red(j) = 1; udca = udupdt(udca,red,df(i)*Ric{i}(j,j)); end % over j lrica=lrica+df(i)*lric(i,:); end % over components end %if end % iterations %recalculating lric into the triangular matrix % for i=1:ncom % ......................... cycle over all components % ric=Ric{i}; [r,d]=ld2ld(udca); pom=lrica'; xx=zeros(nPsi -1,1); if d(nouts+1,nouts+1)>eps, xx(nouts+1)=pom(nouts+1)/d(nouts+1,nouts+1)/2; else xx(nouts+1)=0; end for j=nouts+2:nPsi-1, ff=r(nouts+1:end-1,j)'*d(nouts+1:end-1,nouts+1:end-1)*2*xx(nouts+1:end); if d(j,j)>eps,xx(j)=(pom(j)-ff)/2/d(j,j); else xx(j)=0; end end udca(:,end)= [xx ;0]; % end %puting calculated control factors into a mixture for i=1:ncom, % cl= Ric{i}(nouts+1:nychn,:); for j= nouts+1:nychn, jj=aMix.states.strc(1,j); Fac=facarxls(jj,strc(:,j+1:end)); Fac.Eth=-udca(j,j+1:end); if udca(j,j)