root/applications/dual/SIDP/transformace(3)/sidp_transformace.asv @ 1105

Revision 1105, 4.3 kB (checked in by zimamiro, 14 years ago)
Line 
1function sidp_transformace(sidp_parameters, rsss_parameters,system,apriori)
2%transformace
3N=sidp_parameters.horizont-2;
4poc=zeros(N,2);
5kon=zeros(N,2);
6
7eps=10^-2;
8system.dim=2;
9
10if abs(apriori.y0)-2*apriori.y0_range<0
11    apriori.eta0=(abs(apriori.y0)+2*apriori.y0_range)/2/system.sigma;
12    apriori.eta0_range=apriori.eta0;
13else
14    apriori.eta0=abs(apriori.y0)/system.sigma;
15    apriori.eta0_range=2*apriori.y0_range/system.sigma;
16end
17etar=system.yr/system.sigma;
18
19apriori.beta0abs(apriori.b0)/sqrt(apriori.P0)
202*apriori.b0_range/sqrt(apriori.P0-apriori.P0_range)
21if -<0
22    apriori.beta0=(abs(apriori.b0)+2*apriori.b0_range)/2/sqrt(apriori.P0-apriori.P0_range);
23    apriori.beta0_range=apriori.beta0-eps;
24else
25    apriori.beta0=abs(apriori.b0)/sqrt(apriori.P0);
26    apriori.beta0_range=2*apriori.by0_range/system.sigma;
27end
28
29apriori.beta0=eps;
30apriori.beta0_range=(abs(apriori.b0) + 2*apriori.b0_range)/sqrt(apriori.P0-apriori.P0_range);
31
32%pomocne promenne a konstanty
33step=[1; sidp_parameters.n_grid];
34
35[H ny]=init_hyperstate(N,sidp_parameters.n_grid,apriori,0);  %vytvori H a mi, nastavi H0 a mi*
36
37for i=1:sidp_parameters.n_pass
38   
39    for j=1:sidp_parameters.n_iter
40        [i j]
41        search_region=sidp_parameters.gama^(j-1)*sidp_parameters.lambda^(i-1)*sidp_parameters.search_region_init;
42        [H ny]=update_hyperstate(H, ny);   %vygeneruje trajektorie, v zasazenem regionu rovnomerne rozmisti body a prekopiruje nalezena mi*
43      %>1
44     
45      %los=mc_study2(system,apriori,1000);
46      %porovnej(los)
47     
48        poc(:,:)=H(:,1,:);
49        kon(:,:)=H(:,end,:);     
50        scale=(kon-poc)/(sidp_parameters.n_grid-1);
51     
52        for k=N:-1:1    %pozor na meze
53            level=N-k+1;
54            k
55           
56            for l=1:size(H,2) 
57                candidates=generate_candidates(ny(k,l), search_region, sidp_parameters.num_of_candidates,sidp_parameters.generate_candidates_mode);
58             
59                %compare candidates
60                realization1=zeros(size(candidates,1),rsss_parameters.n0);
61               
62                for m=1:sidp_parameters.num_of_candidates
63                    best_control=-candidates(m);
64                    for n=1:rsss_parameters.n0
65                        eta=H(k,l,1);
66                        beta=H(k,l,2);
67                        %generate realization
68                        for o=1:level
69                            s=randn;
70                            pom=sqrt(1+best_control^2);
71                            eta=abs(eta+beta*best_control+pom*s);
72                            beta=abs(pom*beta+best_control*s);
73                           
74                            realization1(m,n)=realization1(m,n)+(eta-etar(k+o))^2;
75                           
76                            if o<level
77                                %find in H
78                                stav=round(([eta beta]-poc(k+o,:))./scale(k+o,:)+1);
79                                index=min(max(stav,1),sidp_parameters.n_grid)-[0 1];
80                                index=index*step;
81                               
82                                best_control=-ny(k+o,index);
83                                %index==find_in_hyperstate5(state, H(k+o,:,:))
84                                %plot(H(k+o,:,1),H(k+o,:,2), 'r.' ); hold on; plot(eta,beta, 'b+' ); plot(H(k+o,index,1),H(k+o,index,2), 'g+' ); hold off
85                            end
86                        end
87                        best_control=-(eta-etar(end))*beta/(1+beta^2);
88                        s=randn;
89                        pom=sqrt(1+best_control^2);
90                        eta=abs(eta+beta*best_control+pom*s);
91                       
92                        realization1(m,n)=realization1(m,n)+(eta-etar(end))^2;
93                    end
94                end
95                    mean_values=mean(realization1,2);
96                    [min_val min_index]=min(mean_values);
97                    best_control=candidates(min_index);
98                    ny(k,l)=best_control;
99            end
100        end
101       
102        vypis(:,:)=H(:,:,1);
103        save 'eta.txt' vypis -ASCII;
104        vypis(:,:)=H(:,:,2);
105        save 'Beta.txt' vypis -ASCII;
106        save 'ny.txt' ny -ASCII;
107       
108       %  los=mc_study2(system,apriori,1000);
109       %  porovnej(los(1,:))
110    end
111end
112
113end
Note: See TracBrowser for help on using the browser.