ControllerRobust/Main.jl at master · Roberock/ControllerRobust · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
using Pkg
import Pkg;
Pkg.add("StatsPlots")
Pkg.add("ControlSystems");
Pkg.add("Distributed");
Pkg.add("StatsBase");
Pkg.add("Distributions");
Pkg.add("CSV")
# Pkg.add("ControlToolbox")


using ControlSystems, StatsPlots
using Distributed, StatsBase, Distributions
using CSV, DataFrames

include(string(pwd(), "\\src\\build_plant.jl"));
include(string(pwd(), "\\src\\g_controller.jl"));
include(string(pwd(), "\\src\\tf2ss.jl"));
include(string(pwd(), "\\src\\MonteCarloController.jl"));
include(string(pwd(), "\\src\\EstimateReliabilityScores.jl"));
include(string(pwd(), "\\data\\MyDGM.jl"));## my Data-Generating-MechanismDGM
#####################################################################################
####################  Define contorl problem and run example ########################
#####################################################################################

## Parameter δ (uncertain parameters)
δref=[1, 1, 1, 0, 1,0.001]      ## example: a reference value for the uncertain parameters δ
δnames=["mass1" "mass2" "kl" "kn" "λ" "τ" ]     ## names of the uncertain factors
# admissible ranges
δlims=[0.1    2;        0.1    2;
       0.05   1.75;     -1     1;
       0.2    1.8;      0.0001 0.3];

# nominal design
dnom=[-0.1324, 0.3533, 0.6005, 0.0728, 0.5503, 1.4175, 2.6531, 2.4802, 1];
d_names =["b3" "b2" "b1" "b0" "a3" "a2" "a1" "a0"]     ## names of design parameters

### Example: build the system matrices (state-space representation)
A, B, C, D = build_plant(δref);

### Example: define the controller convering from trasfer function to system state matrices
Ac, Bc, Cc, Dc=tf2ss(dnom[1:4],dnom[5:end]);

### Example: call solver and compute reliability function and system response
tvec=0:0.01:25;
g1_nom,g2_nom,g3_nom, y_nominal =g_controller(dnom,δref,tvec,true);

## Example: Monte Carlo simulation  (for random samples within δlims)
Ns=1500
δuniform =  δlims[:,1] .+ (δlims[:,2]- δlims[:,1]).*rand(6,Ns);
# δgaussian=  MyDGM(δref, δlims, Ns)

Gref ,Yref =  MonteCarloController(dnom,Array(δuniform'),tvec)

 ## Exanoke plot system response
plot(tvec,  Yref' , alpha=0.5 ,legend = false, xlabel="time", ylabel="x2(t)")
plot(tvec[1:500],  Yref[:,1:500]' , alpha=0.5 ,legend = false, xlabel="time", ylabel="x2(t)")

## Example:  reliability scores from given-data
Pf_ind, Pf_all, Sev_ind, Sev_all, W, IsF_all, IsF_ind= EstimateReliabilityScores(Gref);


######################################################################
####################  Analyze one of 4 designs #######################
######################################################################
# Candidate controller design
#   𝒮𝒫1(𝒟,0) 𝒮𝒫1(𝒟,.05) 𝒮𝒫2(𝒟) Nominal
# 𝑎4  0.2238   0.5375    0.7600  0.5503
# 𝑎3  0.6811   1.3346    1.9491  1.4175
# 𝑎2  3.1275   2.4206    3.0497  2.6531
# 𝑎1  2.3615   2.1689    2.7344  2.4802
# 𝑎0  1.1833   0.8084    1.0594  1.0000
# 𝑏3 −0.0982   2.4802   −0.0831 -0.1324
# 𝑏2  0.4702   0.6146    0.6358  0.3533
# 𝑏1  0.5886   0.5265    0.7752  0.6005
# 𝑏0  0.0777   0.0716    0.0981  0.0728

Design_Selection=1
if Design_Selection==1 # same as the nominal
    d_new=[-0.1324, 0.3533, 0.6005, 0.0728, 0.5503, 1.4175, 2.6531, 2.4802, 1];
elseif Design_Selection==2 # desing from min{ \hat{F}_w^{-1}(1) }
    d_new=[ -0.103539257684432,	0.452027845175231,	0.553521292381580,	0.0731453358298691,	0.285291345952459, 0.703095280659925,	3.45687805905316,	2.49234017067302,	1.16716573966267]; # 0.717520994282477
elseif Design_Selection==3 # from min { \hat{F}_w^{-1}(0.95) }
    d_new=[-0.0276450664037319, 0.613233915821064, 0.528302902131644, 0.0732150560642133, 0.585865170685746, 1.35042413219186, 2.51075569540054, 2.17663965704309, 0.827018013978000];
elseif Design_Selection==4 # from  min { Pf }
    d_new=[-0.0831455337335365,	0.635784535253079,	0.775222237017883,	0.0981477664191349,	0.760033487719863,	1.94907797910329,	3.04966829435795,	2.73439891631291,	1.05935063251108]; #  gamma: 0.2218
end

## load data (samples of the uncertain factors)

DataSetType =1
Nuber_of_samples=10^3
if DataSetType==1
    File_path=string(pwd(), "\\data\\Theta_small.csv"); ## small uncertainty
    df = CSV.read(File_path,DataFrame);
    δ_samples=Matrix(df);
elseif DataSetType==2
    File_path=string(pwd(), "\\data\\Theta_large.csv"); ## large uncertainty
    df = CSV.read(File_path,DataFrame);
    δ_samples=Matrix(df);
#elseif DataSetType==3
#    δ_samples=MyDGM(Nuber_of_samples); % your own data generating mechanism
end


## Monte Carlo simulation for d_new, failure probability, and severity
G, Y_samples  =  MonteCarloController(d_new,δ_samples,tvec)

Pf_ind, Pf_all, Sev_ind, Sev_all, W, IsFail_all, IsF_ind= EstimateReliabilityScores(G);

Y_safe=Y_samples[vec(IsFail_all.==0),:]';
Y_fail=Y_samples[vec(IsFail_all),:]';
 ## plot system response
plot(tvec,  Y_safe , alpha=0.1 ,legend = false, xlabel="time", ylabel="x2(t)")
plot!(tvec, Y_fail , color =[:red], legend = false, xlabel="time", ylabel="x2(t)")

## plot samples of g1,g2,g3
 l = @layout [a ; b c]

 p1 = scatter( G[vec(IsFail_all),1], G[vec(IsFail_all),2],labels="fail" ,markercolor=[:red],xlabel="g1", ylabel="g2")
 p1 = scatter!( G[vec(IsFail_all.==0),1], G[vec(IsFail_all.==0),2],labels="safe" ,markercolor=[:blue],xlabel="g1", ylabel="g2")

 p2 = scatter( G[vec(IsFail_all),1], G[vec(IsFail_all),3],labels="fail" ,markercolor=[:red],xlabel="g1", ylabel="g3")
 p2 = scatter!( G[vec(IsFail_all.==0),1], G[vec(IsFail_all.==0),3],labels="safe" ,markercolor=[:blue],xlabel="g1", ylabel="g3")

 p3 = scatter( G[vec(IsFail_all),2], G[vec(IsFail_all),3],labels="fail" ,markercolor=[:red],xlabel="g2", ylabel="g3")
 p3 = scatter!( G[vec(IsFail_all.==0),2], G[vec(IsFail_all.==0),3],labels="safe" ,markercolor=[:blue],xlabel="g2", ylabel="g3")

 plot(p1, p2, p3, layout = l)