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

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