[567] | 1 | function [aMix] = soptim(aMix, aMixu, ufc, nstep, chis)
|
---|
| 2 | % soptim performs simultaneous advisory design for normal mixture
|
---|
| 3 | %
|
---|
| 4 | % [aMix] = soptim(aMix, aMixu, ufc, nstep, chis)
|
---|
| 5 | % [aMix] = soptim(aMix, aMixu, ufc, nstep) chis = 1
|
---|
| 6 | % [aMix] = soptim(aMix, aMixu, ufc) nstep = [200, 1]
|
---|
| 7 | %
|
---|
| 8 | % aMix : advised mixture of the type ARX LS enriched on following control states:
|
---|
| 9 | % strc : common control structure
|
---|
| 10 | % ufc : normalised vector qualifying components:
|
---|
| 11 | % dangerous component (0), not dangerous (positive number)
|
---|
| 12 | % kc : lift of quadratic forms
|
---|
| 13 | % UDc : cell vector of u'du decompositions of KLD kernels
|
---|
| 14 | % udca : u'du decomposition of average KLD kernel in UDc
|
---|
| 15 | % kca : average lift of quadratic forms kc
|
---|
| 16 | % aMixu : desired mixture (user's target) of the type ARX LS with control states
|
---|
| 17 | % ufc : vector qualifying components: 0 - dangerous component, (1) - not
|
---|
| 18 | % nstep : parameters [ns1,per] determining design horizon, i.e. horizon = ns1*per;
|
---|
| 19 | % ns1 : number of block repetition
|
---|
| 20 | % per : horizon of a block
|
---|
| 21 | % if nstep is defined by parameter nsl only then per is set to 1
|
---|
| 22 | % chis : indicates strategy chosen: chis=1 for receding horizon (default) and chis=-1 for IST
|
---|
| 23 | %
|
---|
| 24 | % Design : J. Bohm
|
---|
| 25 | % Updated : June, 2002
|
---|
| 26 | % Project : ProDaCTools, IST-1999-12058
|
---|
| 27 | % See also : udupdt, getdvect, facchng, facarxls
|
---|
| 28 |
|
---|
| 29 | % References : \ref{ch9}
|
---|
| 30 | % Note :
|
---|
| 31 | % Updated :
|
---|
| 32 |
|
---|
| 33 |
|
---|
| 34 | 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
|
---|
| 35 |
|
---|
| 36 | if (nargin <4), nstep=[200,1];chis=1;end
|
---|
| 37 | if (nargin <5), chis=1;end
|
---|
| 38 |
|
---|
| 39 | if(isempty(aMix.states.uchn)) error('uoptim needs nonempty list of channels with recognisable actions'); end;
|
---|
| 40 | % normalisation of ufc
|
---|
| 41 | ufc = ufc/sum(ufc);
|
---|
| 42 |
|
---|
| 43 | %Inititialization
|
---|
| 44 | ncom = length(aMix.dfcs);
|
---|
| 45 | dfcs = aMix.dfcs;
|
---|
| 46 | strc = aMix.states.strc; % common control structure
|
---|
| 47 | nPsi = max(size(strc)); % length of regression vector + data
|
---|
| 48 | pochn = aMix.states.pochn; % list of channels with o-innovations
|
---|
| 49 | npochn = length(pochn); % number of channels with o-innovations
|
---|
| 50 | nychn = length(aMix.states.modelled); % number of modelled channels
|
---|
| 51 | nouts = length(aMix.states.outs); % number of innovation channels
|
---|
| 52 | npsi = nPsi-nychn; % length of the regression vector
|
---|
| 53 | kc0 = aMixu.states.kc;
|
---|
| 54 | udca = aMixu.states.udca;
|
---|
| 55 | kca = aMixu.states.kca;
|
---|
| 56 | coms = aMix.coms;
|
---|
| 57 | Ethz = zeros(1,nPsi);
|
---|
| 58 | UDc = aMixu.states.UDc;
|
---|
| 59 | lss = length(nstep);
|
---|
| 60 | lrica = zeros(1,nPsi-1);
|
---|
| 61 | %test of aMixu
|
---|
| 62 | coves=zeros(1,npochn);
|
---|
| 63 | % if npochn~=length(aMixu.Facs),error('aMixu not correctly set'); end
|
---|
| 64 | for i=1:npochn,
|
---|
| 65 | coves(i)=aMixu.Facs{i}.cove;
|
---|
| 66 | end
|
---|
| 67 | if ~any(coves), error('aMixu not correctly set, Facs{.}.cove must be >0'),end
|
---|
| 68 |
|
---|
| 69 | for i=1:ncom, % cycle over number of components ncom
|
---|
| 70 | Ric{i}= zeros(nPsi); % KLD kernels
|
---|
| 71 | kcc(i)=0; % lift of quadratic forms
|
---|
| 72 | end
|
---|
| 73 | df=dfcs/sum(dfcs);
|
---|
| 74 | % setting of design horizon
|
---|
| 75 | if lss==2
|
---|
| 76 | steps = nstep(1)*nstep(2); per=nstep(2);
|
---|
| 77 | else
|
---|
| 78 | steps = nstep; per=1;
|
---|
| 79 | end
|
---|
| 80 | if chis>0, % if the strategy starts from zero
|
---|
| 81 | udca = zeros(nPsi); kca =0;
|
---|
| 82 | end
|
---|
| 83 | %Main design cycle, iterations over the horizon of the criterion
|
---|
| 84 | for iter=1:steps, % ===================== iterations till design horizon
|
---|
| 85 | % ricmn is an auxiliary array accumulating results of optimization
|
---|
| 86 | if mod(iter-1,per)==0,
|
---|
| 87 | % shift of a matrix from bottom right to top left by nychn
|
---|
| 88 | % if nPsi>nychn+1
|
---|
| 89 | [udca, lrica]= ricshift(udca,lrica,nychn,nPsi,npsi);
|
---|
| 90 | % end %if nPsi
|
---|
| 91 | udca(nPsi,nPsi)= 0;
|
---|
| 92 | end %if mod
|
---|
| 93 |
|
---|
| 94 | for i=1:ncom % ......................... cycle over all components
|
---|
| 95 | if mod(iter-1,per)==0,
|
---|
| 96 | % if ufc(i)==0, kcc(i)=1e30; continue; end % excluding bad components
|
---|
| 97 | ric = udca; % ric is auxiliarry working array
|
---|
| 98 | lric(i,:)=lrica;
|
---|
| 99 | kcc(i) = -npochn +kc0(i);
|
---|
| 100 | % adding to each component its stationary loss
|
---|
| 101 | for j=1:npsi, %
|
---|
| 102 | red=UDc{i}(j,:);
|
---|
| 103 | red(j)=1;
|
---|
| 104 | ric= udupdt(ric,red,ufc(i)*UDc{i}(j,j));
|
---|
| 105 | end % for j
|
---|
| 106 | else
|
---|
| 107 | % iteration continues in corresponding component kernel and lift
|
---|
| 108 | % shift of a matrix from bottom right to top left by nychn
|
---|
| 109 | if nPsi>nychn+1,
|
---|
| 110 | [ric,lric(i,:)]=ricshift(Ric{i},lric(i,:),nychn,nPsi,npsi);
|
---|
| 111 | end % if nPsi
|
---|
| 112 | ric(nPsi,nPsi) = Ric{i}(nPsi,nPsi);
|
---|
| 113 | end %if mod
|
---|
| 114 | % expectation is calculated channel by channel
|
---|
| 115 | for j=1: nouts %--------------------- cycle over innovation channels
|
---|
| 116 | indv = ~isempty(find(strc(1,j)==pochn)); % indicator if the channel is o-innovation
|
---|
| 117 | [ric,lric,kcc]= ricexp(ric,lric,kcc,i,j,aMix,nPsi);
|
---|
| 118 | % the penalization is used
|
---|
| 119 | if indv, % visibility indicator
|
---|
| 120 | [ric,lric,kcc]= ricpen(ric,lric,kcc,i,j,aMix,aMixu,nPsi);
|
---|
| 121 | end %end if indv
|
---|
| 122 | end % ------------------------------ reduced all factors of i-th component
|
---|
| 123 | % now penalization
|
---|
| 124 | for j=nouts+1: nychn,
|
---|
| 125 | [ric,lric,kcc]= ricpenu(ric,lric,kcc,i,j,aMix,aMixu,nPsi);
|
---|
| 126 | % if ric(j,j)>eps,
|
---|
| 127 | % kcc(i)=kcc(i)-lric(i,j)*lric(i,j)/ric(j,j)/4;
|
---|
| 128 | % disp('pred m');
|
---|
| 129 | % keyboard
|
---|
| 130 | % % lric(i,:)=lric(i,:)-lric(i,j)*[zeros(1,j) ric(j,j+1:end-1)];
|
---|
| 131 | % disp('po m');
|
---|
| 132 | % keyboard
|
---|
| 133 | % end % if
|
---|
| 134 | end %end for
|
---|
| 135 |
|
---|
| 136 | [l,d]= ld2ld(ric(nouts+1: nychn,nouts+1: nychn));
|
---|
| 137 | 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);
|
---|
| 138 | Ric{i} = ric;
|
---|
| 139 | end % ...................... done for all components
|
---|
| 140 | % now putting losses together
|
---|
| 141 | if mod(iter,per)==0,
|
---|
| 142 | udca=zeros(size(ric));
|
---|
| 143 | lrica = zeros(1,nPsi-1);
|
---|
| 144 | det=ones(1,ncom);
|
---|
| 145 | %%putting it together
|
---|
| 146 | for i=1:ncom,
|
---|
| 147 | % if ufc(i)==0, continue; end
|
---|
| 148 | for j=nouts+1:nPsi,
|
---|
| 149 | red = Ric{i}(j,:);
|
---|
| 150 | red(j) = 1;
|
---|
| 151 | udca = udupdt(udca,red,df(i)*Ric{i}(j,j));
|
---|
| 152 | end % over j
|
---|
| 153 | lrica=lrica+df(i)*lric(i,:);
|
---|
| 154 |
|
---|
| 155 | end % over components
|
---|
| 156 | end %if
|
---|
| 157 |
|
---|
| 158 | end % iterations
|
---|
| 159 | %recalculating lric into the triangular matrix
|
---|
| 160 | % for i=1:ncom % ......................... cycle over all components
|
---|
| 161 | % ric=Ric{i};
|
---|
| 162 | [r,d]=ld2ld(udca);
|
---|
| 163 | pom=lrica';
|
---|
| 164 | xx=zeros(nPsi -1,1);
|
---|
| 165 | if d(nouts+1,nouts+1)>eps, xx(nouts+1)=pom(nouts+1)/d(nouts+1,nouts+1)/2;
|
---|
| 166 | else xx(nouts+1)=0;
|
---|
| 167 | end
|
---|
| 168 | for j=nouts+2:nPsi-1,
|
---|
| 169 | ff=r(nouts+1:end-1,j)'*d(nouts+1:end-1,nouts+1:end-1)*2*xx(nouts+1:end);
|
---|
| 170 | if d(j,j)>eps,xx(j)=(pom(j)-ff)/2/d(j,j);
|
---|
| 171 | else xx(j)=0;
|
---|
| 172 | end
|
---|
| 173 | end
|
---|
| 174 | udca(:,end)= [xx ;0];
|
---|
| 175 | % end
|
---|
| 176 |
|
---|
| 177 | %puting calculated control factors into a mixture
|
---|
| 178 | for i=1:ncom,
|
---|
| 179 | % cl= Ric{i}(nouts+1:nychn,:);
|
---|
| 180 | for j= nouts+1:nychn,
|
---|
| 181 | jj=aMix.states.strc(1,j);
|
---|
| 182 | Fac=facarxls(jj,strc(:,j+1:end));
|
---|
| 183 | Fac.Eth=-udca(j,j+1:end);
|
---|
| 184 | if udca(j,j)<eps, error('synthesis uoptim uncorrect, check its input parameters'),end
|
---|
| 185 | Fac.cove=1/udca(j,j);
|
---|
| 186 | aMix=facchng(aMix,i,Fac);
|
---|
| 187 | end
|
---|
| 188 | end
|
---|
| 189 | aMix.states.udca = udca;
|
---|
| 190 | aMix.states.kca = kca;
|
---|
| 191 | aMix.states.UDc = Ric;
|
---|
| 192 | aMix.states.kc = kcc;
|
---|
| 193 |
|
---|