root/applications/bdmtoolbox/tutorial/mpdm/dist_ctrl_acc.m @ 883

Revision 883, 2.2 kB (checked in by smidl, 14 years ago)

bdmtoolbox + arx

Line 
1clear all;
2% name random variables
3y1 = RV({'y1'},1);
4y2 = RV({'y2'},1);
5y3 = RV({'y3'},1);
6u1 = RV({'u1'},1);
7u2 = RV({'u2'},1);
8
9% create f(y_t| y_{t-3}, u_{t-1})
10fy.class = 'mlnorm<ldmat>';
11fy.rv    = RVjoin([y1,y2,y3]);
12fy.rvc   = RVtimes([y1,y2,y3,u1,u1,u2,u2], [-1, -1, -1, 0, -1, 0, -1]);
13fy.A     = [0.8 , 0.2 , 0 , -0.3 , 0.4 , 0 , 0;...
14           -0.2 , 0.5 , -0.8 , 0.2 , 0.5 , -0.2 , -0.5;...
15            0 , 1.1 , -0.5 , 0 , 0 , -0.2 , 0.3];
16fy.const = [0;0;0];
17fy.R     = 0.1*eye(3);
18
19DS.class = 'PdfDS';
20DS.pdf = fy;
21
22% create ARX estimator
23A1.class = 'ARX';
24A1.rv = RVjoin([y1,y2]);
25A1.rgr = RVtimes([y1,y2,u1,u1],[-1, -1, 0, -1]) ; % correct structure is {y,y}
26A1.log_level ='logbounds,logevidence';
27A1.frg = 0.99;
28
29A2=A1;
30A2.rv = RVjoin([y2,y3]);
31A2.rgr = RVtimes([y2,y3,u2,u2],[-1, -1, 0, -1]) ; % correct structure is {y,y}
32
33C1.class = 'LQG_ARX';
34C1.ARX = A1;
35C1.Qu = 0.1;
36C1.Qy = 0.1*eye(2);
37C1.yreq = [0;1]; %y2=1
38C1.horizon = 1;
39
40C2=C1;
41C2.ARX = A2;
42C2.yreq = [1;0]; %y2=1
43
44P1.class = 'ARXAgent';
45P1.name = 'P1';
46P1.lqg_arx = C1;
47P1.lqg_arx.class = 'LQG_ARX';
48P1.merger.class = 'merger_mix';
49P1.merger.method = 'geometric';
50%P1.merger.dbg_file = 'mp.it';
51P1.merger.ncoms = 1;
52P1.merger.stop_niter= 5;
53P1.neighbours = {};%{'P2'};
54
55P2=P1;
56P2.name = 'P2';
57P2.lqg_arx = C2;
58P2.neighbours = {};
59
60exper.Ndat = 10;
61exper.burnin = 3;
62exper.burn_pdf.class = 'enorm<ldmat>';
63exper.burn_pdf.mu = [0;0];
64exper.burn_pdf.R  = 0.01*eye(2);
65
66
67%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MONTE CARLO %%%%%%%%%%%%%%%%%%%
68
69Ntrials = 100;
70loss_non_coop = zeros(1,Ntrials);
71for i=1:Ntrials
72    M= arena(DS,{P1,P2},exper);
73
74    Y = [M.DS_y1 M.DS_y2 M.DS_y3];
75    Yreq = ones(size(M.DS_y1))*[0 1 0];
76    loss_non_coop(i) = trace((Y-Yreq)'*0.01*(Y-Yreq)) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
77    if loss_non_coop(i)>100
78        %keyboard
79    end
80end
81mean(loss_non_coop)
82
83loss_coop = zeros(1,Ntrials);
84for i=1:Ntrials
85    P1.neighbours = {'P2'};
86    P2.neighbours = {'P1'};
87    M= arena(DS,{P1,P2},exper);
88
89    Y = [M.DS_y1 M.DS_y2 M.DS_y3];
90    Yreq = ones(size(M.DS_y1))*[0 1 0];
91    loss_coop(i) = trace((Y-Yreq)'*0.01*(Y-Yreq)) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
92    if loss_coop(i)>100
93        %keyboard
94    end
95end
96mean(loss_coop)
97
Note: See TracBrowser for help on using the browser.