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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
function varargout = compare(varargin)
[lhs,rhs] = argn(0)
//------------------------------------------------------------------------------
// checking the number of inputs
if rhs < 2 || rhs > 3 then
error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"compare"))
end
//------------------------------------------------------------------------------
data = varargin(1)
model = varargin(2)
if rhs == 3 then
kStep = varargin(3)
elseif rhs == 2 then
kStep = %inf
end
//------------------------------------------------------------------------------
// k step analysis
if typeof(kStep) <> 'constant' || isnan(kStep) then
error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"compare"))
end
// if given k step is infinity or []
// if isinf(kStep) || ~size(kStep,'*') then
// kStep = 1
// end
// checking the dimensions
// if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then
// error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict"))
// end
//------------------------------------------------------------------------------
// checking the plant model
if typeof(model) ~= 'idpoly' then
error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"compare"))
end
modelSampleTime = model.Ts
modelTimeUnit = model.TimeUnit
//------------------------------------------------------------------------------
//checking the data type
if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then
error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"compare"))
end
// checking the plant data
if typeof(data) == 'iddata' then
if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then
error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare"))
end
plantSampleTime = data.Ts
plantTimeUnit = data.TimeUnit
data = [data.OutputData data.InputData]
//disp('iddata')
elseif typeof(data) == 'constant' then
if size(data,'c') ~= 2 then
error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare"))
end
plantSampleTime = model.Ts
plantTimeUnit = model.TimeUnit
end
//------------------------------------------------------------------------------
// comparing the sampling time
if modelSampleTime-plantSampleTime <> 0 then
error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"compare"))
end
// Comparing the time units
if ~strcmp(modelTimeUnit,plantTimeUnit) then
else
error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"compare"))
end
//------------------------------------------------------------------------------
// ckecking the k step size. if it greater than number of sample size then the
// k step will become 1
// if kStep >= size(data,'r') then
// kStep = 1
// end
//------------------------------------------------------------------------------
//storing the plant data
// B(z) C(z)
// y(n) = ---------- u(n) + ---------- e(n)
// A(z)*F(z) A(z)*D(z)
aPoly = poly(model.a,'q','coeff');
bPoly = poly(model.b,'q','coeff');
cPoly = poly(model.c,'q','coeff');
dPoly = poly(model.d,'q','coeff');
fPoly = poly(model.f,'q','coeff');
Gq = bPoly/(aPoly*fPoly)
Hq = cPoly/(aPoly*dPoly)
//disp(kStep)
if kStep == %inf then
//disp('in inf')
outData = sim(data(:,2),model)
else
if kStep == 1 then
Wkq = Hq^-1
elseif kStep > 1 then
adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1);
adPoly = poly(adCoeff,'q','coeff')
cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1);
cPoly = poly(cCoeff,'q','coeff')
hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001)
hBarPoly = poly(hBar,'q','coeff')
Wkq = hBarPoly*Hq^-1//*bPoly/(dPoly*fPoly)
end
// pause
WkqGq = Wkq * Gq
tempWkqGq = coeff(WkqGq.den)
if tempWkqGq(1) <> 1 then
WkqGq.num = WkqGq.num/tempWkqGq(1)
WkqGq.den = WkqGq.den/tempWkqGq(1)
end
Wkq1 = 1-Wkq
tempWkq1 = coeff(Wkq1.den)
if tempWkq1(1) == 1 then
Wkq1.num = Wkq1.num/tempWkq1(1)
Wkq1.den = Wkq1.den/tempWkq1(1)
end
z = poly(0,'z');
WkqGqNum = horner(WkqGq.num,1/z);WkqGqDen = horner(WkqGq.den,1/z);
uPoly = WkqGqNum/WkqGqDen;
uPoly = syslin('d',uPoly);uData = flts(data(:,2)',uPoly);
Wkq1Num = horner(Wkq1.num,1/z);Wkq1Den = horner(Wkq1.den,1/z);
yPoly = Wkq1Num/Wkq1Den;
yPoly = syslin('d',yPoly);yData = flts(data(:,1)',yPoly);
outData = (uData+yData)'
end
tData = (1:size(data,'r'))'*plantSampleTime
fitData = fitValue(data(:,1),outData)
if lhs == 1 then
varargout(1) = []
plot(tData,data(:,1),'m')
plot(tData,outData,'b')
legend('Plant Data','Model Data : '+string(fitData)+'%')
xtitle('Comparison of Time Response','Time ('+ plantTimeUnit+')','Amplitude')
xgrid()
elseif lhs == 2 then
varargout(1) = fitData
varargout(2) = outData
elseif lhs == 3 then
varargout(1) = outData
varargout(2) = tData
varargout(3) = fitData
end
endfunction
|
