|
| 1 | +% Author: Mehran Attar |
| 2 | +% Written: 08-March-2023 |
| 3 | +% Last update: |
| 4 | +% Last revision:--- |
| 5 | + |
| 6 | +%------------- BEGIN CODE -------------- |
| 7 | +clc |
| 8 | +clear all |
| 9 | +close all |
| 10 | +%% system dynamics |
| 11 | +A=[0.7969 -0.2247; |
| 12 | + 0.1798 0.9767]; |
| 13 | +B=[0.1271;0.0132]; |
| 14 | +C = [1,0]; |
| 15 | +D = 0; |
| 16 | +rand('seed',1); |
| 17 | +% define continuous time system |
| 18 | +% sys_c = ss(A,B,C,D); |
| 19 | +% % convert to discrete system |
| 20 | +% samplingtime = 0.01; |
| 21 | +% sys_d = c2d(sys_c,samplingtime); |
| 22 | +dim_x = size(A, 2); |
| 23 | +initpoints =1; |
| 24 | +%Number of time steps |
| 25 | +steps = 4; |
| 26 | +totalsamples = initpoints*steps; |
| 27 | +%% initial set and input |
| 28 | +% X0 = zonotope(zeros(dim_x,1),diag(ones(dim_x,1))); |
| 29 | +% X0 = zonotope(zeros(2,1),0.106*[0.999 -0.1;-1 -1]); |
| 30 | +U = zonotope(0,3); |
| 31 | +X0 = zonotope(zeros(2,1),10*eye(2)); |
| 32 | +%noise zontope W |
| 33 | +W = zonotope(zeros(dim_x,1),0.005*eye(dim_x)); |
| 34 | + |
| 35 | +%Construct matrix zonotpe \mathcal{M}_w |
| 36 | +index=1; |
| 37 | +for i=1:size(W.generators,2) |
| 38 | + vec = W.Z(:,i+1); |
| 39 | + GW{index}= [vec,zeros(dim_x,totalsamples-1)]; |
| 40 | + for j=1:totalsamples-1 |
| 41 | + GW{j+index}= [GW{index+j-1}(:,2:end) GW{index+j-1}(:,1)]; |
| 42 | + end |
| 43 | + index = j+index+1; |
| 44 | +end |
| 45 | +Wmatzono= matZonotope(zeros(dim_x,totalsamples),GW); |
| 46 | + |
| 47 | +% randomly choose constant inputs for each step / sampling time |
| 48 | +for i=1:totalsamples |
| 49 | + u(i) = randPoint(U); |
| 50 | +end |
| 51 | + |
| 52 | +%simulate the system to get the data |
| 53 | +x0 = X0.center; |
| 54 | +x(:,1) = x0; |
| 55 | +index=1; |
| 56 | +for j=1:dim_x:initpoints*dim_x |
| 57 | + x(j:j+dim_x-1,1) = randPoint(X0); |
| 58 | + for i=1:steps |
| 59 | + utraj(j,i) = u(index); |
| 60 | + x(j:j+dim_x-1,i+1) = A*x(j:j+dim_x-1,i) + B*u(index) + randPoint(W); |
| 61 | + index=index+1; |
| 62 | + end |
| 63 | +end |
| 64 | + |
| 65 | +% concatenate the data trajectories |
| 66 | +index_0 =1; |
| 67 | +index_1 =1; |
| 68 | +for j=1:dim_x:initpoints*dim_x |
| 69 | + for i=2:steps+1 |
| 70 | + x_meas_vec_1(:,index_1) = x(j:j+dim_x-1,i); |
| 71 | + index_1 = index_1 +1; |
| 72 | + end |
| 73 | + for i=1:steps |
| 74 | + u_mean_vec_0(:,index_0) = utraj(j,i); |
| 75 | + x_meas_vec_0(:,index_0) = x(j:j+dim_x-1,i); |
| 76 | + index_0 = index_0 +1; |
| 77 | + end |
| 78 | +end |
| 79 | + |
| 80 | +% X_+ is X_1T |
| 81 | +% X_- is X_0T |
| 82 | +U_full = u_mean_vec_0(:,1:totalsamples); %same as u |
| 83 | +X_0T = x_meas_vec_0(:,1:totalsamples); |
| 84 | +X_1T = x_meas_vec_1(:,1:totalsamples); |
| 85 | + |
| 86 | +X1W_cen = X_1T - Wmatzono.center; |
| 87 | +X1W = matZonotope(X1W_cen,Wmatzono.generator); |
| 88 | +% set of A and B |
| 89 | +AB = X1W * pinv([X_0T;U_full]); |
| 90 | + |
| 91 | +%% |
| 92 | + |
| 93 | +matrixCenter = [AB.center;zeros(1,3)]; |
| 94 | + |
| 95 | +for i=1:AB.gens |
| 96 | + G{i}=[AB.generator{i};zeros(1,3)]; |
| 97 | +end |
| 98 | + |
| 99 | +% instantiate matrix zonotope |
| 100 | +M_zono = matZonotope(matrixCenter, G); |
| 101 | + |
| 102 | +% obtain result of all vertices------------------------- |
| 103 | +V_AB = vertices(M_zono); |
| 104 | + |
| 105 | +save V_AB |
0 commit comments