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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
//[y,x0] = predict(data,idpoly,k)
//References
//Digital Control(11.1.2) by Kanna M.Moudgalya
//System Identification Theory for User Second Edition (3.2) by Lennart Ljung
// Code Aurthor - Ashutosh Kumar Bhargava
function varargout = predict(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"),"predict"))
end
//------------------------------------------------------------------------------
data = varargin(1)
model = varargin(2)
if rhs == 3 then
kStep = varargin(3)
elseif rhs == 2 then
kStep = 1
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"),"predict"))
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"),"predict"))
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"),"predict"))
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"),"predict"))
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"),"predict"))
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"),"predict"))
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"),"predict"))
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)
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
end
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
//pause
//------------------------------------------------------------------------------
// storing the plant data
uCoeff = coeff(WkqGq.num*Wkq1.den)
yCoeff = coeff(WkqGq.den*Wkq1.num)
yCapCoeff = coeff(WkqGq.den*Wkq1.den)
//pause
lengthuCoeff = length(uCoeff)
lengthyCoeff = length(yCoeff)
lengthyCapCoeff = length(yCapCoeff)
//------------------------------------------------------------------------------
// keeping initial conditions equal to zero
uData = zeros(lengthuCoeff,1)
yData = zeros(lengthyCoeff,1)
yCapData = zeros(lengthyCapCoeff-1,1)
uData = [uData;data(:,2)]
yData = [yData;data(:,1)]
sampleData = size(data,'r')
//pause
// reversing the coefficients
if ~size(uCoeff,'*') then
uCoeff = 0
else
uCoeff = uCoeff(lengthuCoeff:-1:1)
end
if ~size(yCoeff) then
yCoeff = 0
else
yCoeff = yCoeff(lengthyCoeff:-1:1)
end
if ~size(yCapCoeff,'*') then
yCapCoeff = 0
else
yCapCoeff = -yCapCoeff(lengthyCapCoeff:-1:2)
end
//pause
for ii = 1:sampleData+1
//pause
if ~size(uData(ii:ii+lengthuCoeff-1),'*') then
tempu = 0
else
tempu = uCoeff*uData(ii:ii+lengthuCoeff-1);
end
if ~size(yData(ii:ii+lengthyCoeff-1),'*')
tempy = 0
else
tempy = yCoeff*yData(ii:ii+lengthyCoeff-1);
end
if ~size(yCapData(ii:ii+lengthyCapCoeff-2),'*') then
tempyCap = 0
else
tempyCap = yCapCoeff*yCapData(ii:ii+lengthyCapCoeff-2);
end
yCapData = [yCapData;tempu+tempy+tempyCap];
end
// pause
extraSample = abs(size(yCapData,'r')-sampleData)
yCapData = yCapData(extraSample+1:$)
timeData = ((modelSampleTime:modelSampleTime:(size(yCapData,'r')*modelSampleTime))');
if lhs == 1 then
clf()
plot(timeData,yCapData)
axisData = gca()
tempTimeUnit = 'Time('+modelTimeUnit+')'
xtitle('Predicted Response',tempTimeUnit,'y')
xgrid
varargout(1) = 0
elseif lhs == 2 then
varargout(1) = yCapData
varargout(2) = 0
elseif lhs == 3 then
varargout(1) = yCapData
varargout(2) = timeData
varargout(3) = 0
end
endfunction
|
