Changeset 755
- Timestamp:
- 12/14/09 20:36:06 (14 years ago)
- Location:
- applications/dual/IterativeLocal
- Files:
-
- 1 added
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
applications/dual/IterativeLocal/ildp.m
r748 r755 10 10 % [b=0.6, yr=10] sigma = 0.1; rho = 1.0; 11 11 12 Iterace = 20; %iterace13 K = 20; %casy12 Iterace = 10; %iterace 13 K = 5; %casy 14 14 N = 100; %vzorky 15 15 … … 36 36 37 37 Xkn = zeros(3, K, N); 38 Xstripe = zeros(3, K); 39 40 41 % Kpi = [[1.0308 0.9990 0.9797 0.9899 1.0075 0.9830 1.0050 1.0173 0.9659 1.0000]; 42 % [0.4104 0.5562 0.6933 0.4674 0.3817 0.5036 0.4203 0.3461 0.7750 1.0000]; 43 % [0.9204 0.8378 0.6707 0.3629 0.9815 0.8955 1.2227 1.1550 1.9374 1.0000]; 44 % [0.5756 0.4225 0.4055 0.5941 0.5188 0.3218 0.3435 0.3210 0.0333 0]]; 45 38 Xstripe = zeros(3, K); 39 46 40 gka = 0; 47 41 gnu = 0; … … 49 43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50 44 %iteracni smycka 51 % clf reset 52 53 % PtMin = 0.0001; 54 55 % UU=zeros(K,N); 45 56 46 for i = 1:Iterace, 57 47 %generovani stavu 58 % hold off59 48 for n = 1:N, 60 49 Xkn(:, 1, n) = x0 + [sigma*randn(); 0; 0]; 61 50 for k = 1:K-1, 62 51 Uk = uPi(k, Xkn(:, k, n)); 63 % UU(k,n) = Uk;64 52 Kk = Uk*Xkn(3, k, n)/(Uk^2*Xkn(3, k, n) + sigma^2); 65 53 Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn(); 66 54 Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk); 67 55 Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n); 68 % plot(1:K,Xkn(1,:,n)) 69 end 70 % hold all 71 % subplot(4,1,1); 72 % plot(1:K,Xkn(1,:,n)) 73 % hold all 74 % subplot(4,1,2); 75 % plot(1:K,Xkn(2,:,n)) 76 % hold all 77 % subplot(4,1,3); 78 % plot(1:K,Xkn(3,:,n)) 79 % hold all 80 % subplot(4,1,4); 81 % plot(1:K-1,UU(:,n)) 82 56 end 83 57 end 84 58 Xstripe = mean(Xkn, 3); 85 % hold all 86 % subplot(4,1,1); 87 % plot(1:K,Xstripe(1,:),'-ro') 88 % hold all 89 % subplot(4,1,2); 90 % plot(1:K,Xstripe(2,:),'-ro') 91 % hold all 92 % subplot(4,1,3); 93 % plot(1:K,Xstripe(3,:),'-ro') 94 95 %hold off 96 % Epsl = randn(3,N); 97 %iterace pro k = K-1..1 59 98 60 for k = K-1:-1:1, 99 61 gka = k; … … 103 65 Xkn(1, k, n) = Xstripe(1, k) + rho*randn(); 104 66 Xkn(2, k, n) = Xstripe(2, k) + rho*randn(); 105 Xkn(3, k, n) = Xstripe(3, k)*exp(rho*randn()); 106 % %rovne okoli 107 % Xkn(1, k, n) = Xstripe(1, k) + Epsl(1,n); 108 % Xkn(2, k, n) = Xstripe(2, k) + Epsl(2,n); 109 % Xkn(3, k, n) = Xstripe(3, k)*exp(Epsl(3,n)); 110 end 111 % Xkn(:,k,:) 112 67 Xkn(3, k, n) = Xstripe(3, k)*exp(rho*randn()); 68 end 113 69 % 2] 114 70 for n = 1:N, 115 71 gnu = n; 116 % if(k == 1)117 % Uopt(n) = Rk(k)/2;118 % Hmin(n) = Hamilt(Uopt(n));119 % else120 72 [Uopt(n), Hmin(n)] = fminunc(@Hamilt, uPi(k, Xkn(:, k, n)), optimset('GradObj','on','Display','notify')); 121 % % [i, k, n]122 % end123 %124 % interv = -1000:1:1000;125 % for ll = 1:2001,126 % hodnot(ll) = Hamilt(interv(ll));127 % end128 % plot(interv, hodnot,'-b',Uopt(n),Hmin(n),'rs',uPi(k, Xkn(:, k, n)),Hamilt(uPi(k, Xkn(:, k, n))),'go')129 % prah = 100;130 % if(Uopt(n) > prah)131 % Uopt(n) = prah;132 % Hmin(n) = Hamilt(prah);133 % disp('u > horni mez');134 % elseif(Uopt(n) < -prah)135 % Uopt(n) = -prah;136 % Hmin(n) = Hamilt(-prah);137 % disp('u < dolni mez');138 % end139 % if(extfl < 1)140 % disp('exitflag < 1')141 % end142 73 end 143 74 % 3] 144 for n = 1:N, % V??? nema to byt k+1? ale asi ne 145 Vn(n) = Hmin(n) + Vtilde(k+1, Xkn(:, k, n)); 146 % xxx(n) = Xkn(1,k,n); 147 % xxy(n) = Xkn(2,k,n); 148 % xxz(n) = Xkn(3,k,n); 149 % Vn(n) = xxx(n).*xxx(n).*xxx(n); 150 end 151 % xxx2 = sort(xxx); 152 % xxy2 = sort(xxy); 153 % xxz2 = sort(xxz); 154 75 for n = 1:N, 76 Vn(n) = Hmin(n) + Vtilde(k+1, Xkn(:, k, n)); 77 end 155 78 % 4] 156 79 Epsilon = zeros(3, N); … … 163 86 FiFiTInvFi = (mFi*mFi')\mFi; 164 87 Wv(:,k) = FiFiTInvFi * Vn'; 165 % Wv = zeros(10,1);166 % Wv(1, k) = Wvtmp(1);167 % Wv(2, k) = Wvtmp(2);168 % Wv(3, k) = Wvtmp(3);169 % Wv(4, k) = Wvtmp(4);170 % FiFiTInvFi = (mFi'*mFi)\mFi';171 % Wv(:, k) = FiFiTInvFi' * Vn';172 % UU(k,:) = Uopt;173 % for n = 1:N,174 % rozd(n) = Vn(n) - Vtilde(k,Xkn(:,k,n));175 % end176 177 88 for n = 1:N, 178 89 yt(n) = Xkn(1, k, n); … … 185 96 Uopt']; 186 97 PsiPsiTInvPsi = (mPsi'*mPsi)\mPsi'; 187 Kpi(:,k) = PsiPsiTInvPsi * (Rk(k)*ones(N,1)); 188 189 % for nn=1:N, 190 % vlv(nn) = Vtilde(k, [xxx2(nn);xxy2(nn);xxz2(nn)]); 191 % end 192 % subplot(3,1,1); 193 % plot(xxx,Vn,'rs',xxx2,vlv,'-b') 194 % subplot(3,1,2); 195 % plot(xxy,Vn,'rs',xxy2,vlv,'-b') 196 % subplot(3,1,3); 197 % plot(xxz,Vn,'rs',xxz2,vlv,'-b') 198 % clf reset 199 % hold off 200 % 201 % yii = 0.5:0.025:1.5; 202 % bjj = 0.5:0.025:1.5; 203 % for ii= 1:41, 204 % for jj= 1:41, 205 % vmtrx(ii,jj) = Vtilde(k, [yii(ii);bjj(jj);0]); 206 % end 207 % end 208 % 209 % % xlabel('yt') 210 % % ylabel('bt') 211 % % zlabel('Vt') 212 % surf(yii,bjj,vmtrx) 213 % % % for n=1:N, 214 % % % hold all 215 % % % surf(yt(n),bt(n),Vn(n),'LineStyle','','Marker','o'); 216 % % % end 217 % hold all 218 % plot3(bt, yt, Vn, 'ro') 219 220 98 Kpi(:,k) = PsiPsiTInvPsi * (Rk(k)*ones(N,1)); 221 99 end 222 % clf reset223 % for n=1:N,224 %225 % hold all226 % subplot(4,1,1);227 % plot(1:K,Xkn(1,:,n),'-b')228 % hold all229 % subplot(4,1,2);230 % plot(1:K,Xkn(2,:,n),'-b')231 % hold all232 % subplot(4,1,3);233 % plot(1:K,Xkn(3,:,n),'-b')234 % hold all235 % % subplot(4,1,4);236 % % plot(1:K,UU(:,n),'-b')237 % end238 % for k=1:K,239 % riz(k) = uPi(k, Xstripe(:, k));240 % ce(k) = (Rk(k) - Xstripe(1, k))/(Xstripe(2, k) + Xstripe(3, k));241 % end242 % plot(1:K,riz,1:K,ce,1:K,Xstripe(1,:))243 244 100 end 245 101 %%%%%%%%%%% … … 247 103 248 104 Kpi 249 %graficky vystup 250 251 % X1 = zeros(3, K); 252 % UU1 = zeros(1,K); 253 % 254 % X1(:,1) = x0 + [sigma*randn(); 0; 0]; 255 % for k = 1:K-1, 256 % Upi = uPi(k, X1); 257 % UU1(k) = Upi; 258 % Ktmp = Upi*X1(3,k)/(Upi^2*X1(3,k) + sigma^2); 259 % X1(1,k+1) = X1(1,k)+X1(2,k)*Upi + sigma*randn(); 260 % X1(2,k+1) = X1(2,k) + Ktmp*(X1(1,k+1) - X1(1,k) - X1(2,k)*Upi); 261 % X1(3,k+1) = X1(3,k)*(1-Ktmp*Upi); 262 % end 263 % % X 264 % % hold off 265 % subplot(4,1,1); 266 % plot(1:K,X(1,:),'-gs') 267 % % hold off 268 % subplot(4,1,2); 269 % plot(1:K,X(2,:),'-gs') 270 % % hold off 271 % subplot(4,1,3); 272 % plot(1:K,X(3,:),'-gs') 273 % % hold off 274 % subplot(4,1,4); 275 % plot(1:K,UU,'-gs') 276 % 277 % figure 278 % for k=1:K, 279 % riz(k) = uPi(k, Xstripe(:, k)); 280 % ce(k) = (Rk(k) - Xstripe(1, k))/(Xstripe(2, k) + Xstripe(3, k)); 281 % end 282 % plot(1:K,riz,1:K,ce,1:K,Xstripe(1,:)) 105 283 106 for n = 1:N, 284 107 Xkn(:, 1, n) = x0 + [sigma*randn(); 0; 0]; … … 289 112 Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn(); 290 113 Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk); 291 Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n); 292 % plot(1:K,Xkn(1,:,n)) 114 Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n); 293 115 end 294 116 hold all … … 316 138 Xkn(1, k+1, n) = Xkn(1, k, n) + Xkn(2, k, n)*Uk + sigma*randn(); 317 139 Xkn(2, k+1, n) = Xkn(2, k, n) + Kk*(Xkn(1, k+1, n) - Xkn(1, k, n) - Xkn(2, k, n)*Uk); 318 Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n); 319 % plot(1:K,Xkn(1,:,n)) 140 Xkn(3, k+1, n) = (1 - Kk*Uk)*Xkn(3, k, n); 320 141 end 321 142 hold all … … 347 168 hist(log(loss)) 348 169 349 % disp()350 % X = zeros(1, K);351 % b = 0;352 % UU = zeros(1,K);353 %354 % X(1) = sigma*randn();355 % for k = 1:K-1,356 % Upi = uPi(k, [X,b,0]);357 % UU(k) = Upi;358 % X(k+1) = X(k) + b*Upi + sigma*randn();359 % end360 %361 % % hold off362 % subplot(2,1,1);363 % plot(1:K,X,'-gs')364 % % hold off365 % subplot(2,1,2);366 % plot(1:K,UU,'-gs')367 368 170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 369 171 %pomocne funkce … … 416 218 end 417 219 418 % function [val_Vt] = Vtilde_dx_dx(k_Vt, x_Vt)419 % if(k_Vt == K)420 % val_Vt = hdxdx;421 % else422 % val_Vt = zeros(3,3);423 % val_Vt(1,1) = 2*Wv(5,k_Vt);424 % val_Vt(2,2) = 2*Wv(8,k_Vt);425 % val_Vt(3,3) = 2*Wv(10,k_Vt);426 %427 % val_Vt(1,2) = Wv(6,k_Vt);428 % val_Vt(1,3) = Wv(7,k_Vt);429 % val_Vt(2,3) = Wv(9,k_Vt);430 %431 % val_Vt(2,1) = val_Vt(1,2);432 % val_Vt(3,1) = val_Vt(1,3);433 % val_Vt(3,2) = val_Vt(2,3);434 % end435 % end436 437 % function [val_Vt] = trSgVt(k_Vt)438 % if(k_Vt == K)439 % val_Vt = 0;440 % else441 % val_Vt = 2*Wv(5,k_Vt)*sigma;442 % end443 % end444 445 220 function [val_Fi] = vectFi(x_Fi) 446 221 val_Fi = [ ... … … 469 244 x_Fi(2, :).^2; ... 470 245 x_Fi(2, :).*log(x_Fi(3, :)); ... 471 % 2*ln(x_Fi(3, :)); ...472 246 ]; 473 % val_Fi = [ ... 474 % ones(1, N); ... 475 % x_Fi(1, :); ... 476 % x_Fi(2, :); ... 477 % x_Fi(3, :); ... 478 % ]; 247 479 248 end 480 249 … … 490 259 0 2*x_Fi(2) 0; ... 491 260 0 log(x_Fi(3)) x_Fi(2)/(x_Fi(3)); ... 492 % 0 0 2*x_Fi(3); ...493 261 ]; 494 262 end
