source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/control/BlackBox.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 94.7 KB
Line 
1/* Class BlackBox
2*
3* This class contains the constructor to create an instance of
4* a generalised BlackBox with a single input, single output
5* and a gain. It contins the methods for obtaining the
6* transfer function in the s-domain and the z-domain.
7*
8* This class is the superclass for several sub-classes,
9* e.g. Prop (P controller), PropDeriv (PD controller),
10* PropInt (PI controller), PropIntDeriv (PID controller),
11* FirstOrder, SecondOrder, AtoD (ADC), DtoA (DAC),
12* ZeroOrderHold, DelayLine, OpenLoop (Open Loop Path),
13* of use in control engineering.
14*
15* Author: Michael Thomas Flanagan.
16*
17* Created: August 2002
18* Updated: 17 July 2003, 18 May 2005, 6 April 2008, 6 October 2009, 30 October 2009,
19* 2-9 November 2009, 20 January 2010, 23-25 May 2010, 3 June 2010, 18 January 2011
20*
21*
22* DOCUMENTATION:
23* See Michael T Flanagan's JAVA library on-line web page:
24* http://www.ee.ucl.ac.uk/~mflanaga/java/BlackBox.html
25* http://www.ee.ucl.ac.uk/~mflanaga/java/
26*
27* Copyright (c) 2002 - 2011 Michael Thomas Flanagan
28*
29* PERMISSION TO COPY:
30*
31* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
32* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
33* and associated documentation or publications.
34*
35* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, this list of conditions
36* and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
37*
38* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
39* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
40*
41* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
42* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
43* or its derivatives.
44*
45***************************************************************************************/
46
47
48package agents.anac.y2015.agentBuyogV2.flanagan.control;
49
50import agents.anac.y2015.agentBuyogV2.flanagan.complex.*;
51import agents.anac.y2015.agentBuyogV2.flanagan.interpolation.CubicSpline;
52import agents.anac.y2015.agentBuyogV2.flanagan.io.Db;
53import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
54import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
55import agents.anac.y2015.agentBuyogV2.flanagan.plot.Plot;
56import agents.anac.y2015.agentBuyogV2.flanagan.plot.PlotGraph;
57import agents.anac.y2015.agentBuyogV2.flanagan.plot.PlotPoleZero;
58
59
60public class BlackBox{
61
62 protected int sampLen = 0; // Length of array of stored inputs, outputs and time
63 protected double[] inputT = null; // Array of input signal in the time domain
64 protected double[] outputT = null; // Array of output signal in the time domain
65 protected double[] time = null; // Array of time at which inputs were taken (seconds)
66 protected double forgetFactor = 1.0D; // Forgetting factor, e.g. in exponential forgetting of error values
67 protected double deltaT = 0.0D; // Sampling time (seconds)
68 protected double sampFreq = 0.0D; // Sampling frequency (Hz)
69 protected Complex inputS = new Complex(); // Input signal in the s-domain
70 protected Complex outputS = new Complex(); // Output signal in the s-domain
71 protected Complex sValue = new Complex(); // Laplacian s
72 protected Complex zValue = new Complex(); // z-transform z
73 protected ComplexPoly sNumer = new ComplexPoly(1.0D); // Transfer function numerator in the s-domain
74 protected ComplexPoly sDenom = new ComplexPoly(1.0D); // Transfer function denominator in the s-domain
75 protected ComplexPoly zNumer = new ComplexPoly(1.0D); // Transfer function numerator in the z-domain
76 protected ComplexPoly zDenom = new ComplexPoly(1.0D); // Transfer function denominator in the z-domain
77 protected boolean sNumerSet = false; // = true when numerator entered
78 protected boolean sDenomSet = false; // = true when denominator entered
79 protected Complex sNumerScaleFactor = Complex.plusOne();// s-domain numerator/(product of s - zeros)
80 protected Complex sDenomScaleFactor = Complex.plusOne();// s-domain denominator/(product of s - poles)
81 protected Complex sNumerWorkingFactor = Complex.plusOne();// s-domain numerator/(product of s - zeros) at that point in the program
82 protected Complex sDenomWorkingFactor = Complex.plusOne();// s-domain denominator/(product of s - poles) at that point in the program
83 protected Complex[] sPoles = null; // Poles in the s-domain
84 protected Complex[] sZeros = null; // Zeros in the s-domain
85 protected Complex[] zPoles = null; // Poles in the z-domain
86 protected Complex[] zZeros = null; // Zeros in the z-domain
87 protected int sNumerDeg = 0; // Degree of transfer function numerator in the s-domain
88 protected int sDenomDeg = 0; // Degree of transfer function denominator in the s-domain
89 protected int zNumerDeg = 0; // Degree of transfer function numerator in the z-domain
90 protected int zDenomDeg = 0; // Degree of transfer function denominator in the z-domain
91 protected double deadTime = 0.0D; // Time delay between an input and the matching output [in s-domain = exp(-s.deadTime)]
92 protected int orderPade = 2; // Order(1 to 4)of the pade approximation for exp(-sT)
93 // default option = 2
94 protected ComplexPoly sNumerPade = new ComplexPoly(1.0D); // Transfer function numerator in the s-domain including Pade approximation
95 protected ComplexPoly sDenomPade = new ComplexPoly(1.0D); // Transfer function denominator in the s-domain including Pade approximation
96 protected Complex[] sPolesPade = null; // Poles in the s-domain including Pade approximation
97 protected Complex[] sZerosPade = null; // Zeros in the s-domain including Pade approximation
98 protected int sNumerDegPade = 0; // Degree of transfer function numerator in the s-domain including Pade approximation
99 protected int sDenomDegPade = 0; // Degree of transfer function denominator in the s-domain including Pade approximation
100 protected boolean maptozero = true; // if true infinity s zeros map to zero
101 // if false infinity s zeros map to minus one
102 protected boolean padeAdded = false; // if true Pade poles and zeros added
103 // if false No Pade poles and zeros added
104 protected double integrationSum=0.0D; // Stored integration sum in numerical integrations
105 protected int integMethod = 1; // numerical integration method
106 // = 0 Trapezium Rule [default option]
107 // = 1 Backward Rectangular Rule
108 // = 2 Foreward Rectangular Rule
109 protected int ztransMethod = 0; // z trasform method
110 // = 0 s -> z mapping (ad hoc procedure) from the continuous time domain erived s domain functions
111 // = 1 specific z transform, e.g. of a difference equation
112 protected String name = "BlackBox"; // Superclass or subclass name, e.g. pid, pd, firstorder.
113 // user may rename an instance of the superclass or subclass
114 protected String fixedName = "BlackBox"; // Super class or subclass permanent name, e.g. pid, pd, firstorder.
115 // user must NOT change fixedName in any instance of the superclass or subclass
116 // fixedName is used as an identifier in classes such as OpenPath, ClosedLoop
117 protected int nPlotPoints = 400; // number of points used tp lot response curves, e.g. step input response curve
118
119 protected String[] subclassName = {"BlackBox", "OpenLoop", "ClosedLoop", "Prop", "PropDeriv", "PropInt", "PropIntDeriv", "FirstOrder", "SecondOrder", "Compensator", "LowPassPassive", "HighPassPassive", "Transducer", "DelayLine", "ZeroOrderHold", "AtoD", "DtoA"};
120 protected int nSubclasses = subclassName.length; // number of subclasses plus superclass
121 protected int subclassIndex = 0; // = 0 BlackBox
122 // = 1 OpenLoop
123 // = 2 ClosedLoop
124 // = 3 Prop
125 // = 4 PropDeriv
126 // = 5 PropInt
127 // = 6 PropIntDeriv
128 // = 7 FirstOrder
129 // = 8 SecondOrder
130 // = 9 Compensator
131 // = 10 LowPassPassive
132 // = 11 HighPassPassive
133 // = 12 Transducer
134 // = 13 DelayLine
135 // = 14 ZeroOrderHold
136 // = 15 AtoD
137 // = 16 DtoA
138
139 // Constructor
140 public BlackBox(){
141 }
142
143 // Constructor with fixedName supplied
144 // for use by subclasses
145 public BlackBox(String name){
146 this.name = name;
147 this.fixedName = name;
148 this.setSubclassIndex();
149 }
150
151 // Set subclass index
152 protected void setSubclassIndex(){
153 boolean test = true;
154 int i = 0;
155 while(test){
156 if(this.fixedName.equals(subclassName[i])){
157 this.subclassIndex = i;
158 test = false;
159 }
160 else{
161 i++;
162 if(i>=this.nSubclasses){
163 System.out.println("Subclass name, " + this.fixedName + ", not recognised as a recorder subclass");
164 System.out.println("Subclass, " + this.fixedName + ", handled as BlackBox");
165 this.subclassIndex = i;
166 test = false;
167 }
168 }
169 }
170 }
171
172 // Set the transfer function numerator in the s-domain
173 // Enter as an array of real (double) coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
174 public void setSnumer(double[] coeff){
175 this.sNumerDeg = coeff.length-1;
176 this.sNumer = new ComplexPoly(coeff);
177 this.sNumerSet = true;
178 this.calcPolesZerosS();
179 this.addDeadTimeExtras();
180 }
181
182 // Method to set extra terms to s-domain numerator and denominator and
183 // to calculate extra zeros and poles if the dead time is not zero.
184 protected void addDeadTimeExtras()
185 {
186 this.sNumerDegPade = this.sNumerDeg;
187 this.sNumerPade = this.sNumer.copy();
188 this.sDenomDegPade = this.sDenomDeg;
189 this.sDenomPade = this.sDenom.copy();
190
191 if(this.deadTime==0.0D){
192 this.transferPolesZeros();
193 }
194 else{
195 this.pade();
196 }
197
198 }
199
200 // Set the transfer function numerator in the s-domain
201 // Enter as an array of Complex coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
202 public void setSnumer(Complex[] coeff){
203 this.sNumerDeg = coeff.length-1;
204 this.sNumer = new ComplexPoly(coeff);
205 this.sNumerSet = true;
206 this.calcPolesZerosS();
207 this.addDeadTimeExtras();
208
209
210 }
211
212 // Set the transfer function numerator in the s-domain
213 // Enter as an existing instance of ComplexPoly
214 public void setSnumer(ComplexPoly coeff){
215 this.sNumerDeg = coeff.getDeg();
216 this.sNumer = ComplexPoly.copy(coeff);
217 this.sNumerSet = true;
218 this.calcPolesZerosS();
219 this.addDeadTimeExtras();
220 }
221
222 // Set the transfer function denominator in the s-domain
223 // Enter as an array of real (double) coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
224 public void setSdenom(double[] coeff){
225 this.sDenomDeg = coeff.length-1;
226 this.sDenom = new ComplexPoly(coeff);
227 this.sDenomSet = true;
228 this.calcPolesZerosS();
229 this.addDeadTimeExtras();
230 }
231
232 // Set the transfer function denomonator in the s-domain
233 // Enter as an array of Complex coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
234 public void setSdenom(Complex[] coeff){
235 this.sDenomDeg = coeff.length-1;
236 this.sDenom = new ComplexPoly(coeff);
237 this.sDenomSet = true;
238 this.calcPolesZerosS();
239 this.addDeadTimeExtras();
240
241 }
242
243 // Set the transfer function denominator in the s-domain
244 // Enter as an existing instance of ComplexPoly
245 public void setSdenom(ComplexPoly coeff){
246 this.sDenomDeg = coeff.getDeg();
247 this.sDenom = coeff.copy();
248 this.sDenomSet = true;
249 this.calcPolesZerosS();
250 this.addDeadTimeExtras();
251 }
252
253 // calculate constant converting product of root terms to the value of the polynomial
254 public static Complex scaleFactor(ComplexPoly poly, Complex[] roots){
255 int nRoots = roots.length;
256
257 // calculate mean of the poles
258 Complex mean = new Complex(0.0D, 0.0);
259 for(int i=0; i<nRoots; i++)mean = mean.plus(roots[i]);
260 mean = mean.over(nRoots);
261 // check that mean != a root; increase mean by 1.5 till != any pole
262 boolean test = true;
263 int ii=0;
264 while(test){
265 if(mean.isEqual(roots[ii])){
266 if(mean.isEqual(Complex.zero())){
267 for(int i=0; i<nRoots; i++)mean = mean.plus(roots[i].abs());
268 if(mean.isEqual(Complex.zero()))mean=Complex.plusOne();
269 }
270 else{
271 mean = mean.times(1.5D);
272 }
273 ii=0;
274 }
275 else{
276 ii++;
277 if(ii>nRoots-1)test = false;
278 }
279 }
280
281 // calculate product of roots-mean
282 Complex product = new Complex(1.0D, 0.0);
283 for(int i=0; i<nRoots; i++)product = product.times(mean.minus(roots[i]));
284
285 // evaluate the polynomial at mean value
286 Complex eval = poly.evaluate(mean);
287
288 // Calculate scaleFactor
289 return eval.over(product);
290 }
291
292 // Get numerator scale factor
293 public Complex getSnumerScaleFactor(){
294 if(this.sNumerScaleFactor==null)this.calcPolesZerosS();
295 return this.sNumerScaleFactor;
296 }
297
298 // Get denominator scale factor
299 public Complex getSdenomScaleFactor(){
300 if(this.sDenomScaleFactor==null)this.calcPolesZerosS();
301 return this.sDenomScaleFactor;
302 }
303
304 // Set the dead time
305 public void setDeadTime(double deadtime){
306 this.deadTime = deadtime;
307 this.pade();
308 }
309
310 // Set the dead time and the Pade approximation order
311 public void setDeadTime(double deadtime, int orderPade){
312 this.deadTime = deadtime;
313 if(orderPade>5){
314 orderPade=4;
315 System.out.println("BlackBox does not support Pade approximations above an order of 4");
316 System.out.println("The order has been set to 4");
317 }
318 if(orderPade<1){
319 orderPade=1;
320 System.out.println("Pade approximation order was less than 1");
321 System.out.println("The order has been set to 1");
322 }
323 this.orderPade = orderPade;
324 this.pade();
325 }
326
327 // Set the Pade approximation order
328 public void setPadeOrder(int orderPade){
329 if(orderPade>5){
330 orderPade=4;
331 System.out.println("BlackBox does not support Pade approximations above an order of 4");
332 System.out.println("The order has been set to 4");
333 }
334 if(orderPade<1){
335 orderPade=2;
336 System.out.println("Pade approximation order was less than 1");
337 System.out.println("The order has been set to 2");
338 }
339 this.orderPade = orderPade;
340 this.pade();
341 }
342
343 // Get the dead time
344 public double getDeadTime(){
345 return this.deadTime;
346 }
347
348 // Get the Pade approximation order
349 public int getPadeOrder(){
350 return this.orderPade;
351 }
352
353 // Resets the s-domain Pade inclusive numerator and denominator adding a Pade approximation
354 // Also calculates and stores additional zeros and poles arising from the Pade approximation
355 protected void pade(){
356 ComplexPoly sNumerExtra = null;
357 ComplexPoly sDenomExtra = null;
358 Complex[] newZeros = null;
359 Complex[] newPoles = null;
360 switch(orderPade){
361 case 1: this.sNumerDegPade = this.sNumerDeg + 1;
362 this.sDenomDegPade = this.sDenomDeg + 1;
363 this.sNumerPade = new ComplexPoly(sNumerDegPade);
364 this.sDenomPade = new ComplexPoly(sDenomDegPade);
365 sNumerExtra = new ComplexPoly(1.0D, -this.deadTime/2.0D);
366 sDenomExtra = new ComplexPoly(1.0D, this.deadTime/2.0D);
367 this.sNumerPade = this.sNumer.times(sNumerExtra);
368 this.sDenomPade = this.sDenom.times(sDenomExtra);
369 newZeros = Complex.oneDarray(1);
370 newZeros[0].reset(2.0/this.deadTime, 0.0D);
371 newPoles = Complex.oneDarray(1);
372 newPoles[0].reset(-2.0/this.deadTime, 0.0D);
373 break;
374 case 2: this.sNumerDegPade = this.sNumerDeg + 2;
375 this.sDenomDegPade = this.sDenomDeg + 2;
376 this.sNumerPade = new ComplexPoly(sNumerDegPade);
377 this.sDenomPade = new ComplexPoly(sDenomDegPade);
378 sNumerExtra = new ComplexPoly(1.0D, -this.deadTime/2.0D, Math.pow(this.deadTime, 2)/12.0D);
379 sDenomExtra = new ComplexPoly(1.0D, this.deadTime/2.0D, Math.pow(this.deadTime, 2)/12.0D);
380 this.sNumerPade = this.sNumer.times(sNumerExtra);
381 this.sDenomPade = this.sDenom.times(sDenomExtra);
382 newZeros = sNumerExtra.rootsNoMessages();
383 newPoles = sDenomExtra.rootsNoMessages();
384 break;
385 case 3: this.sNumerDegPade = this.sNumerDeg + 3;
386 this.sDenomDegPade = this.sDenomDeg + 3;
387 this.sNumerPade = new ComplexPoly(sNumerDegPade);
388 this.sDenomPade = new ComplexPoly(sDenomDegPade);
389 double[] termn3 = new double[4];
390 termn3[0] = 1.0D;
391 termn3[1] = -this.deadTime/2.0D;
392 termn3[2] = Math.pow(this.deadTime, 2)/10.0D;
393 termn3[3] = -Math.pow(this.deadTime, 3)/120.0D;
394 sNumerExtra = new ComplexPoly(termn3);
395 this.sNumerPade = this.sNumer.times(sNumerExtra);
396 newZeros = sNumerExtra.rootsNoMessages();
397 double[] termd3 = new double[4];
398 termd3[0] = 1.0D;
399 termd3[1] = this.deadTime/2.0D;
400 termd3[2] = Math.pow(this.deadTime, 2)/10.0D;
401 termd3[3] = Math.pow(this.deadTime, 3)/120.0D;
402 sDenomExtra = new ComplexPoly(termd3);
403 this.sDenomPade = this.sDenom.times(sDenomExtra);
404 newPoles = sDenomExtra.rootsNoMessages();
405 break;
406 case 4: this.sNumerDegPade = this.sNumerDeg + 4;
407 this.sDenomDegPade = this.sDenomDeg + 4;
408 this.sNumerPade = new ComplexPoly(sNumerDegPade);
409 this.sDenomPade = new ComplexPoly(sDenomDegPade);
410 double[] termn4 = new double[5];
411 termn4[0] = 1.0D;
412 termn4[1] = -this.deadTime/2.0D;
413 termn4[2] = 3.0D*Math.pow(this.deadTime, 2)/28.0D;
414 termn4[3] = -Math.pow(this.deadTime, 3)/84.0D;
415 termn4[4] = Math.pow(this.deadTime, 4)/1680.0D;
416 sNumerExtra = new ComplexPoly(termn4);
417 this.sNumerPade = this.sNumer.times(sNumerExtra);
418 newZeros = sNumerExtra.rootsNoMessages();
419 double[] termd4 = new double[5];
420 termd4[0] = 1.0D;
421 termd4[1] = this.deadTime/2.0D;
422 termd4[2] = 3.0D*Math.pow(this.deadTime, 2)/28.0D;
423 termd4[3] = Math.pow(this.deadTime, 3)/84.0D;
424 termd4[4] = Math.pow(this.deadTime, 4)/1680.0D;
425 sDenomExtra = new ComplexPoly(termd4);
426 this.sDenomPade = this.sDenom.times(sDenomExtra);
427 newPoles = sDenomExtra.rootsNoMessages();
428 break;
429 default: this.orderPade = 2;
430 this.sNumerDegPade = this.sNumerDeg + 2;
431 this.sDenomDegPade = this.sDenomDeg + 2;
432 this.sNumerPade = new ComplexPoly(sNumerDegPade);
433 this.sDenomPade = new ComplexPoly(sDenomDegPade);
434 sNumerExtra = new ComplexPoly(1.0D, -this.deadTime/2.0D, Math.pow(this.deadTime, 2)/12.0D);
435 sDenomExtra = new ComplexPoly(1.0D, this.deadTime/2.0D, Math.pow(this.deadTime, 2)/12.0D);
436 this.sNumerPade = this.sNumer.times(sNumerExtra);
437 this.sDenomPade = this.sDenom.times(sDenomExtra);
438 newZeros = sNumerExtra.rootsNoMessages();
439 newPoles = sDenomExtra.rootsNoMessages();
440 break;
441 }
442
443 // store zeros and poles arising from the Pade term
444 if(this.sNumerPade!=null && this.sNumerDegPade>0){
445 sZerosPade = Complex.oneDarray(sNumerDegPade);
446 for(int i=0; i<sNumerDeg; i++){
447 sZerosPade[i] = sZeros[i].copy();
448 }
449 for(int i=0; i<this.orderPade; i++){
450 sZerosPade[i+sNumerDeg] = newZeros[i].copy();
451 }
452 }
453
454 if(this.sDenomPade!=null && this.sDenomDegPade>0){
455 sPolesPade = Complex.oneDarray(sDenomDegPade);
456 for(int i=0; i<sDenomDeg; i++){
457 sPolesPade[i] = sPoles[i].copy();
458 }
459 for(int i=0; i<this.orderPade; i++){
460 sPolesPade[i+sDenomDeg] = newPoles[i].copy();
461 }
462 }
463 this.zeroPoleCancellation();
464 this.padeAdded = true;
465 }
466
467 // Copies s-domain poles and zeros from the s-domain arrays to the s-domain Pade arrays
468 // used when deadTime is zero
469 protected void transferPolesZeros(){
470
471 this.sNumerDegPade = this.sNumerDeg;
472 this.sNumerPade = this.sNumer.copy();
473 if(this.sNumerDeg>0 && this.sZeros!=null){
474 this.sZerosPade = Complex.oneDarray(this.sNumerDeg);
475 for(int i=0; i<this.sNumerDeg; i++)this.sZerosPade[i] = this.sZeros[i].copy();
476 }
477
478 this.sDenomDegPade = this.sDenomDeg;
479 this.sDenomPade = this.sDenom.copy();
480 if(this.sDenomDeg>0 && this.sPoles!=null){
481 this.sPolesPade = Complex.oneDarray(this.sDenomDeg);
482 for(int i=0; i<this.sDenomDeg; i++)this.sPolesPade[i] = this.sPoles[i].copy();
483 }
484 this.zeroPoleCancellation();
485 this.padeAdded = true;
486
487 }
488
489 // Get the Pade approximation order
490 public int orderPade(){
491 return this.orderPade;
492 }
493
494 // Warning message if dead time greater than sampling period
495 protected boolean deadTimeWarning(String method){
496 boolean warning = false; // warning true if dead time is greater than the sampling period
497 // false if not
498 if(this.deadTime>this.deltaT){
499 System.out.println(this.name+"."+method+": The dead time is greater than the sampling period");
500 System.out.println("Dead time: "+this.deadTime);
501 System.out.println("Sampling period: "+this.deltaT);
502 System.out.println("!!! The results of this program may not be physically meaningful !!!");
503 warning = true;
504 }
505 return warning;
506 }
507
508 // Perform z transform for a given delta T
509 // Uses maptozAdHoc in this class but may be overridden in a subclass
510 public void zTransform(double deltat){
511 this.mapstozAdHoc(deltat);
512 }
513
514 // Perform z transform using an already set delta T
515 // Uses maptozAdHoc in this class but may be overridden in a subclass
516 public void zTransform(){
517 this.mapstozAdHoc();
518 }
519
520 // Map s-plane zeros and poles of the transfer function onto the z-plane using the ad-hoc method
521 // for a given sampling period.
522 // References:
523 // John Dorsey, Continuous and Discrete Control Systems, pp 490-491, McGraw Hill (2002)
524 // J R Leigh, Applied Digital Control, pp 78-80, Prentice-Hall (1985)
525 public void mapstozAdHoc(double deltaT){
526 this.deltaT = deltaT;
527 this.mapstozAdHoc();
528 }
529
530 // Map s-plane zeros and poles of the transfer function onto the z-plane using the ad-hoc method
531 // for an already set sampling period.
532 // References:
533 // John Dorsey, Continuous and Discrete Control Systems, pp 490-491, McGraw Hill (2002)
534 // J R Leigh, Applied Digital Control, pp 78-80, Prentice-Hall (1985)
535 public void mapstozAdHoc(){
536
537 this.deadTimeWarning("mapstozAdHoc");
538 if(!this.padeAdded)this.transferPolesZeros();
539
540 // Calculate z-poles
541 this.zDenomDeg = this.sDenomDegPade;
542 ComplexPoly root = new ComplexPoly(1);
543 this.zDenom = new ComplexPoly(this.zDenomDeg);
544 if(zDenomDeg>0){
545 this.zPoles = Complex.oneDarray(this.zDenomDeg);
546 for(int i=0; i<this.zDenomDeg; i++){
547 zPoles[i]=Complex.exp(this.sPolesPade[i].times(this.deltaT));
548 }
549 this.zDenom = ComplexPoly.rootsToPoly(zPoles);
550 }
551
552 // Calculate z-zeros
553 // number of zeros from infinity poles
554 int infZeros = this.sDenomDegPade;
555 // check that total zeros does not exceed total poles
556 if(infZeros+this.sNumerDegPade>this.sDenomDegPade)infZeros=this.sDenomDegPade-this.sNumerDegPade;
557 // total number of zeros
558 this.zNumerDeg = this.sNumerDegPade + infZeros;
559 this.zNumer = new ComplexPoly(zNumerDeg);
560 this.zZeros = Complex.oneDarray(zNumerDeg);
561 // zero values
562 if(this.zNumerDeg>0){
563 for(int i=0; i<this.sNumerDegPade; i++){
564 zZeros[i]=Complex.exp(sZerosPade[i].times(this.deltaT));
565 }
566 if(infZeros>0){
567 if(maptozero){
568 for(int i=this.sNumerDegPade; i<this.zNumerDeg; i++){
569 zZeros[i]=Complex.zero();
570 }
571 }
572 else{
573 for(int i=this.sNumerDegPade; i<this.zNumerDeg; i++){
574 zZeros[i]=Complex.minusOne();
575 }
576 }
577 }
578 this.zNumer = ComplexPoly.rootsToPoly(this.zZeros);
579 }
580
581 // Match s and z steady state gains
582 this.sValue=Complex.zero();
583 this.zValue=Complex.plusOne();
584 boolean testzeros = true;
585 while(testzeros){
586 testzeros = false;
587 if(this.sDenomDegPade>0){
588 for(int i=0; i<this.sDenomDegPade; i++){
589 if(this.sPolesPade[i].truncate(3).equals(this.sValue.truncate(3)))testzeros=true;
590 }
591 }
592 if(!testzeros && this.sNumerDegPade>0){
593 for(int i=0; i<this.sDenomDegPade; i++){
594 if(this.sZerosPade[i].truncate(3).equals(this.sValue.truncate(3)))testzeros=true;
595 }
596 }
597 if(!testzeros && this.zDenomDeg>0){
598 for(int i=0; i<this.zDenomDeg; i++){
599 if(this.zPoles[i].truncate(3).equals(this.zValue.truncate(3)))testzeros=true;
600 }
601 }
602 if(!testzeros && this.zNumerDeg>0){
603 for(int i=0; i<this.zDenomDeg; i++){
604 if(this.zZeros[i].truncate(3).equals(this.zValue.truncate(3)))testzeros=true;
605 }
606 }
607 if(testzeros){
608 this.sValue = this.sValue.plus(Complex.plusJay()).truncate(3);
609 this.zValue = Complex.exp(this.sValue.times(this.deltaT).truncate(3));
610 }
611 }
612 Complex gs = this.evalTransFunctS(this.sValue);
613 Complex gz = this.evalTransFunctZ(this.zValue);
614 Complex constant = gs.over(gz);
615 ComplexPoly constantPoly = new ComplexPoly(constant);
616 this.zNumer = this.zNumer.times(constantPoly);
617 }
618
619 // Set the map infinity zeros to zero or -1 option
620 // maptozero: if true infinity s zeros map to zero
621 // if false infinity s zeros map to minus one
622 // default value = false
623 public void setMaptozero(boolean maptozero){
624 this.maptozero = maptozero;
625 }
626
627 // Set the transfer function numerator in the z-domain
628 // Enter as an array of real (double) coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
629 public void setZnumer(double[] coeff){
630 this.zNumerDeg = coeff.length-1;
631 this.zNumer = new ComplexPoly(coeff);
632 this.zZeros = this.zNumer.rootsNoMessages();
633 }
634
635 // Set the transfer function numerator in the z-domain
636 // Enter as an array of Complex coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
637 public void setZnumer(Complex[] coeff){
638 this.zNumerDeg = coeff.length-1;
639 this.zNumer = new ComplexPoly(coeff);
640 this.zZeros = this.zNumer.rootsNoMessages();
641 }
642
643 // Set the transfer function numerator in the z-domain
644 // Enter as an existing instance of ComplexPoly
645 public void setZnumer(ComplexPoly coeff){
646 this.zNumerDeg = coeff.getDeg();
647 this.zNumer = ComplexPoly.copy(coeff);
648 this.zZeros = this.zNumer.rootsNoMessages();
649 }
650
651 // Set the transfer function denominator in the z-domain
652 // Enter as an array of real (double) coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
653 public void setZdenom(double[] coeff){
654 this.zDenomDeg = coeff.length-1;
655 this.zDenom = new ComplexPoly(coeff);
656 this.zPoles = this.zDenom.rootsNoMessages();
657 }
658
659 // Set the transfer function denomonatot in the z-domain
660 // Enter as an array of Complex coefficients of the polynomial a + bs +c.s.s + d.s.s.s + ....
661 public void setZdenom(Complex[] coeff){
662 this.zDenomDeg = coeff.length-1;
663 this.zDenom = new ComplexPoly(coeff);
664 this.zPoles = this.zDenom.rootsNoMessages();
665 }
666
667 // Set the transfer function denominator in the z-domain
668 // Enter as an existing instance of ComplexPoly
669 public void setZdenom(ComplexPoly coeff){
670 this.zDenomDeg = coeff.getDeg();
671 this.zDenom = ComplexPoly.copy(coeff);
672 this.zPoles = this.zDenom.rootsNoMessages();
673 }
674
675 // Set the sampling period
676 public void setDeltaT(double deltaT ){
677 if(this.deltaT==0.0){
678 this.deltaT=deltaT;
679 this.sampFreq=1.0D/this.deltaT;
680 this.deadTimeWarning("setDeltaT");
681 }
682 else{
683
684 String question = "BlackBox setDeltaT: Do you wish to replace the deltaT value, " + this.deltaT + " with " + deltaT;
685 if(Db.yesNo(question)){
686 this.deltaT=deltaT;
687 this.sampFreq=1.0D/this.deltaT;
688 this.deadTimeWarning("setDeltaT");
689 if(this.time!=null){
690 int holdS = this.sampLen;
691 this.sampLen = (int)Math.round(time[this.sampLen-1]/this.deltaT);
692 double[] holdT = Conv.copy(time);
693 double[] holdI = Conv.copy(inputT);
694 this.time = new double[this.sampLen];
695 this.inputT = new double[this.sampLen];
696 CubicSpline cs = new CubicSpline(holdT, holdI);
697 this.time[0] = holdT[0];
698 this.inputT[0] = holdI[0];
699 for(int i=1; i<this.sampLen-1; i++){
700 this.time[i] = this.time[i-1] = this.deltaT;
701 this.inputT[i] = cs.interpolate(this.time[i]);
702 }
703 this.time[sampLen-1] = holdT[holdS];
704 this.inputT[sampLen-1] = holdI[holdS];
705 }
706 }
707 }
708 }
709
710 // Set the forgetting factor
711 public void setForgetFactor(double forget){
712 this.forgetFactor = forget;
713 }
714
715 // Set the sampling frequency
716 public void setSampFreq(double sfreq ){
717 this.sampFreq=sfreq;
718 this.setDeltaT(1.0D/sfreq);
719 }
720
721 // Set the Laplacian s value (s - Complex)
722 public void setS(Complex s){
723 this.sValue = Complex.copy(s);
724 }
725
726 // Set the Laplacian s value (s - real + imaginary parts)
727 public void setS(double sr, double si){
728 this.sValue.reset(sr,si);
729 }
730
731 // Set the Laplacian s value (s - imag, real = 0.0)
732 public void setS(double si){
733 this.sValue.reset(0.0D, si);
734 }
735
736 // Set the z-transform z value (z - Complex)
737 public void setZ(Complex z){
738 this.zValue = Complex.copy(z);
739 }
740
741 // Set the z-transform z value (z - real + imaginary parts)
742 public void setZ(double zr, double zi){
743 this.zValue.reset(zr,zi);
744 }
745
746 // Set the z transform method
747 // 0 = s to z mapping (ad hoc procedure)
748 // 1 = specific z transform, e.g. z transform of a difference equation
749 public void setZtransformMethod(int ztransMethod){
750 if(ztransMethod<0 || ztransMethod>1){
751 System.out.println("z transform method option number " + ztransMethod + " not recognised");
752 System.out.println("z tr methodansform option number set in BlackBox to the default value of 0 (s -> z ad hoc mapping)");
753 this.integMethod = 0;
754 }
755 else{
756 this.ztransMethod = ztransMethod;
757 }
758 }
759
760 // Set the integration method [number option]
761 // 0 = trapezium, 1 = Backward rectangular, 2 = Foreward rectangular
762 public void setIntegrateOption(int integMethod){
763 if(integMethod<0 || integMethod>2){
764 System.out.println("integration method option number " + integMethod + " not recognised");
765 System.out.println("integration method option number set in BlackBox to the default value of 0 (trapezium rule)");
766 this.integMethod = 0;
767 }
768 else{
769 this.integMethod = integMethod;
770 }
771 }
772
773 // Set the integration method [String option]
774 // trapezium; trapezium, tutin. Backward rectangular; back backward. Foreward rectangular; foreward, fore
775 // Continuous time equivalent: continuous, cont
776 public void setIntegrateOption(String integMethodS){
777 if(integMethodS.equals("trapezium") || integMethodS.equals("Trapezium") ||integMethodS.equals("tutin") || integMethodS.equals("Tutin")){
778 this.integMethod = 0;
779 }
780 else{
781 if(integMethodS.equals("backward") || integMethodS.equals("Backward") ||integMethodS.equals("back") || integMethodS.equals("Back")){
782 this.integMethod = 1;
783 }
784 else{
785 if(integMethodS.equals("foreward") || integMethodS.equals("Foreward") ||integMethodS.equals("fore") || integMethodS.equals("Fore")){
786 this.integMethod = 2;
787 }
788 else{
789 System.out.println("integration method option " + integMethodS + " not recognised");
790 System.out.println("integration method option number set in PID to the default value of 0 (trapezium rule)");
791 this.integMethod = 0;
792 }
793 }
794 }
795 }
796
797 // Reset the length of the arrays storing the times, time domain inputs and time domain outputs
798 public void setSampleLength(int samplen){
799 if(samplen==0)throw new IllegalArgumentException("Entered sample length must be greater than zero");
800 if(samplen==1)samplen=2;
801 if(this.sampLen==0){
802 this.sampLen = samplen;
803 this.time = new double[samplen];
804 this.inputT = new double[samplen];
805 this.outputT = new double[samplen];
806 }
807 else{
808 String question = "BlackBox setSampleLength: Do you wish to replace the sample length, " + this.sampLen + " with " + samplen;
809 if(Db.yesNo(question)){
810 int holdS = this.sampLen;
811 this.sampLen=samplen;
812 if(this.time!=null){
813 this.deltaT = this.time[holdS-1]/(samplen-1);
814 double[] holdT = Conv.copy(time);
815 double[] holdI = Conv.copy(inputT);
816 this.time = new double[this.sampLen];
817 this.inputT = new double[this.sampLen];
818 CubicSpline cs = new CubicSpline(holdT, holdI);
819 this.time[0] = holdT[0];
820 this.inputT[0] = holdI[0];
821 for(int i=1; i<this.sampLen-1; i++){
822 this.time[i] = this.time[i-1] = this.deltaT;
823 this.inputT[i] = cs.interpolate(this.time[i]);
824 }
825 this.time[sampLen-1] = holdT[holdS];
826 this.inputT[sampLen-1] = holdI[holdS];
827 }
828 }
829 }
830 }
831
832 // Reset the name of the black box
833 public void setName(String name){
834 this.name=name;
835 }
836
837 // Enter a single current time domain time and input value
838 public void setInputT(double ttime, double inputt){
839 if(this.deltaT==0.0){
840 this.time = new double[2];
841 this.time[0] = 0.0;
842 this.time[1] = ttime;
843 this.inputT = new double[2];
844 this.inputT[0] = inputt;
845 this.inputT[1] = inputt;
846 this.outputT = new double[2];
847 this.sampLen = 2;
848 // this.deltaT = ttime;
849 // this.sampFreq = 1.0/this.deltaT;
850 }
851 else{
852 double delta = this.deltaT;
853 this.sampLen = (int)Math.round(ttime/delta);
854 this.deltaT = ttime/sampLen;
855 if(!Fmath.isEqualWithinLimits(this.deltaT, delta, delta*1e-3)){
856 System.out.println("BlackBox setInputT method; deltaT has been reset from " + delta + " to " + this.deltaT);
857 }
858 this.sampFreq = 1.0/this.deltaT;
859 this.time = new double[this.sampLen];
860 this.time[this.sampLen-1]=ttime;
861 this.inputT = new double[this.sampLen];
862 this.inputT[this.sampLen-1]=inputt;
863 this.outputT = new double[this.sampLen];
864 for(int i=sampLen-2; i>0; i--){
865 this.time[i]= this.time[i+1]-deltaT;
866 this.inputT[i]=inputt;
867 }
868 this.time[0]=0.0;
869 this.inputT[0]=inputt;
870 }
871 }
872
873 // Enter a set of current time domain times and input values
874 public void setInputT(double[] ttime, double[] inputt){
875 int samplen = ttime.length;
876 if(samplen!=inputt.length)throw new IllegalArgumentException("time and input arrays are of different lengths: " + samplen + ", " + inputt.length);
877 if(samplen==1){
878 this.setInputT(ttime[0], inputt[0]);
879 }
880 else{
881 this.sampLen = samplen;
882 this.time = ttime;
883 this.inputT = inputt;
884 this.outputT = new double[this.sampLen];
885 this.deltaT = ttime[this.sampLen]/(this.sampLen - 1);
886 this.sampFreq = 1.0/this.deltaT;
887 }
888 }
889
890
891
892 // Reset s-domain input
893 public void setInputS(Complex input){
894 this.inputS=input;
895 }
896
897 // Reset all inputs, outputs and times to zero
898 public void resetZero(){
899 for(int i=0; i<this.sampLen-1; i++){
900 this.outputT[i] = 0.0D;
901 this.inputT[i] = 0.0D;
902 this.time[i] = 0.0D;
903 }
904 this.outputS = Complex.zero();
905 this.inputS = Complex.zero();
906 this.deltaT = 0.0;
907 this.sampLen = 0;
908 }
909
910 // Calculate the zeros and poles in the s-domain
911 // does not include Pade approximation term
912 protected void calcPolesZerosS(){
913 if(this.sNumer!=null){
914 if(this.sNumer.getDeg()>0)this.sZeros = this.sNumer.rootsNoMessages();
915 if(this.sZeros!=null){
916 this.sNumerScaleFactor = BlackBox.scaleFactor(this.sNumer, this.sZeros);
917 }
918 else{
919 this.sNumerScaleFactor = this.sNumer.coeffCopy(0);
920 }
921 }
922
923 if(this.sDenom!=null){
924 if(this.sDenom.getDeg()>0)this.sPoles = this.sDenom.rootsNoMessages();
925 if(this.sPoles!=null){
926 this.sDenomScaleFactor = BlackBox.scaleFactor(this.sDenom, this.sPoles);
927 }
928 else{
929 this.sDenomScaleFactor = this.sDenom.coeffCopy(0);
930 }
931 }
932 if(this.sNumerPade!=null){
933 if(this.sNumerPade.getDeg()>0)this.sZerosPade = this.sNumerPade.rootsNoMessages();
934 }
935 if(this.sDenomPade!=null){
936 if(this.sDenomPade.getDeg()>0)this.sPolesPade = this.sDenomPade.rootsNoMessages();
937 }
938 }
939
940 // Eliminates identical poles and zeros in the s-domain
941 protected void zeroPoleCancellation(){
942 boolean check = false;
943 boolean testI = true;
944 boolean testJ = true;
945 int i=0;
946 int j=0;
947
948 if(this.sNumerDegPade==0 || this.sDenomDegPade==0)testI=false;
949 if(this.sZerosPade==null || this.sPolesPade==null)testI=false;
950 while(testI){
951 j=0;
952 while(testJ){
953 if(this.sZerosPade[i].isEqual(this.sPolesPade[j])){
954 for(int k=j+1; k<this.sDenomDegPade; k++)this.sPolesPade[k-1] = this.sPolesPade[k].copy();
955 this.sDenomDegPade--;
956 for(int k=i+1; k<this.sNumerDegPade; k++)this.sZerosPade[k-1] = this.sZerosPade[k].copy();
957 this.sNumerDegPade--;
958 check = true;
959 testJ=false;
960 i--;
961 }
962 else{
963 j++;
964 if(j>this.sDenomDegPade-1)testJ=false;
965 }
966 }
967 i++;
968 if(i>this.sNumerDegPade-1)testI=false;
969 }
970 if(check){
971 if(this.sNumerDegPade==0){
972 this.sNumerPade = new ComplexPoly(1.0D);
973 }
974 else{
975 Complex[] holdn = Complex.oneDarray(sNumerDegPade);
976 for(int ii=0; ii<sNumerDegPade; ii++)holdn[i] = this.sZerosPade[ii].copy();
977 this.sZerosPade = holdn;
978 this.sNumerPade = ComplexPoly.rootsToPoly(this.sZerosPade);
979 }
980 if(this.sDenomDegPade==0){
981 this.sDenomPade = new ComplexPoly(1.0D);
982 }
983 else{
984 Complex[] holdd = Complex.oneDarray(sDenomDegPade);
985 for(int ii=0; ii<sDenomDegPade; ii++)holdd[i] = this.sPolesPade[ii].copy();
986 this.sPolesPade = holdd;
987 this.sDenomPade = ComplexPoly.rootsToPoly(this.sPolesPade);
988 }
989 }
990
991 check = false;
992 testI = true;
993 testJ = true;
994 i=0;
995 j=0;
996
997 if(this.sNumerDeg==0 || this.sDenomDeg==0)testI=false;
998 if(this.sZeros==null || this.sPoles==null)testI=false;
999 while(testI){
1000 j=0;
1001 while(testJ){
1002 if(this.sZeros[i].isEqual(this.sPoles[j])){
1003 for(int k=j+1; k<this.sDenomDeg; k++)this.sPoles[k-1] = this.sPoles[k].copy();
1004 this.sDenomDeg--;
1005 for(int k=i+1; k<this.sNumerDeg; k++)this.sZeros[k-1] = this.sZeros[k].copy();
1006 this.sNumerDeg--;
1007 check = true;
1008 testJ=false;
1009 i--;
1010 }
1011 else{
1012 j++;
1013 if(j>this.sDenomDeg-1)testJ=false;
1014 }
1015 }
1016 i++;
1017 if(i>this.sNumerDeg-1)testI=false;
1018 }
1019 if(check){
1020 if(this.sNumerDeg==0){
1021 this.sNumer = new ComplexPoly(1.0D);
1022 }
1023 else{
1024 Complex[] holdn = Complex.oneDarray(sNumerDeg);
1025 for(int ii=0; ii<sNumerDeg; ii++)holdn[i] = this.sZeros[ii].copy();
1026 this.sZeros = holdn;
1027 this.sNumer = ComplexPoly.rootsToPoly(this.sZeros);
1028 this.sNumerWorkingFactor = this.sNumerScaleFactor;
1029 }
1030 if(this.sDenomDeg==0){
1031 this.sDenom = new ComplexPoly(1.0D);
1032 }
1033 else{
1034 Complex[] holdd = Complex.oneDarray(sDenomDeg);
1035 for(int ii=0; ii<sDenomDeg; ii++)holdd[i] = this.sPoles[ii].copy();
1036 this.sPoles = holdd;
1037 this.sDenom = ComplexPoly.rootsToPoly(this.sPoles);
1038 this.sDenomWorkingFactor = this.sDenomScaleFactor;
1039 }
1040 }
1041 }
1042
1043 // Get steadty state value for a unit step input
1044 public double getSeadyStateValue(){
1045 Complex num = this.sNumer.evaluate(Complex.zero());
1046 Complex den = this.sDenom.evaluate(Complex.zero());
1047 Complex ssc = num.over(den);
1048 double ssdr = ssc.getReal();
1049 double ssdi = ssc.getImag();
1050 if(Math.abs(ssdi)>Math.abs(ssdr)*0.01){
1051 System.out.println("method getSteadyStateValue: The imaginary part, " + ssdi + ", is greater than 1 per cent of the the real part, " + ssdr);
1052 System.out.println("Magnitude has been returned");
1053 }
1054 return ssc.abs();
1055 }
1056
1057 // Get steadty state value for a step input of magnitude, mag
1058 public double getSeadyStateValue(double mag){
1059 Complex num = this.sNumer.evaluate(Complex.zero());
1060 Complex den = this.sDenom.evaluate(Complex.zero());
1061 Complex ssc = num.over(den);
1062 double ssdr = ssc.getReal();
1063 double ssdi = ssc.getImag();
1064 if(Math.abs(ssdi)>Math.abs(ssdr)*0.01){
1065 System.out.println("method getSteadyStateValue: The imaginary part, " + ssdi + ", is greater than 1 per cent of the the real part, " + ssdr);
1066 System.out.println("Magnitude has been returned");
1067 }
1068 return mag*ssc.abs();
1069 }
1070
1071
1072
1073 // Evaluate the s-domain tranfer function for the present value of s
1074 // deadtime evaluated as exponential term
1075 public Complex evalTransFunctS(){
1076 if(!this.padeAdded)this.transferPolesZeros();
1077 Complex num = this.sNumer.evaluate(this.sValue);
1078 Complex den = this.sDenom.evaluate(this.sValue);
1079 Complex lagterm = Complex.plusOne();
1080 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1081 return num.over(den).times(lagterm);
1082 }
1083
1084 // Evaluate the s-domain tranfer function for a given Complex value of s
1085 public Complex evalTransFunctS(Complex sValue){
1086 if(!this.padeAdded)this.transferPolesZeros();
1087 this.sValue = Complex.copy(sValue);
1088 Complex num = this.sNumer.evaluate(sValue);
1089 Complex den = this.sDenom.evaluate(sValue);
1090 Complex lagterm = Complex.plusOne();
1091 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1092 return num.over(den).times(lagterm);
1093 }
1094
1095 // Evaluate the s-domain tranfer function for a sine wave input at a given frequency (s^-1)
1096 public Complex evalTransFunctS(double freq){
1097 if(!this.padeAdded)this.transferPolesZeros();
1098 this.sValue.reset(0.0D, 2.0D*Math.PI*freq);
1099 Complex num = this.sNumer.evaluate(this.sValue);
1100 Complex den = this.sDenom.evaluate(this.sValue);
1101 Complex lagterm = Complex.plusOne();
1102 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1103 return num.over(den).times(lagterm);
1104 }
1105
1106 // Evaluate the magnitude of the s-domain tranfer function for the present value of s
1107 public double evalMagTransFunctS(){
1108 if(!this.padeAdded)this.transferPolesZeros();
1109 Complex num = this.sNumer.evaluate(this.sValue);
1110 Complex den = this.sDenom.evaluate(this.sValue);
1111 Complex lagterm = Complex.plusOne();
1112 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1113 return (num.over(den).times(lagterm)).abs();
1114 }
1115
1116 // Evaluate the magnitude of the s-domain tranfer function for a given Complex value of s
1117 public double evalMagTransFunctS(Complex sValue){
1118 if(!this.padeAdded)this.transferPolesZeros();
1119 this.sValue = Complex.copy(sValue);
1120 Complex num = this.sNumer.evaluate(sValue);
1121 Complex den = this.sDenom.evaluate(sValue);
1122 Complex lagterm = Complex.plusOne();
1123 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1124 return (num.over(den).times(lagterm)).abs();
1125 }
1126
1127 // Evaluate the magnitude of the s-domain tranfer function for a sine wave input at a given frequency (s^-1)
1128 public double evalMagTransFunctS(double freq){
1129 if(!this.padeAdded)this.transferPolesZeros();
1130 this.sValue.reset(0.0D, 2.0D*Math.PI*freq);
1131 Complex num = this.sNumer.evaluate(this.sValue);
1132 Complex den = this.sDenom.evaluate(this.sValue);
1133 Complex lagterm = Complex.plusOne();
1134 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1135 return (num.over(den).times(lagterm)).abs();
1136 }
1137
1138 // Evaluate the phase of the s-domain tranfer function for the present value of s
1139 public double evalPhaseTransFunctS(){
1140 if(!this.padeAdded)this.transferPolesZeros();
1141 Complex num = this.sNumer.evaluate(this.sValue);
1142 Complex den = this.sDenom.evaluate(this.sValue);
1143 Complex lagterm = Complex.plusOne();
1144 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1145 return (num.over(den).times(lagterm)).arg();
1146 }
1147
1148 // Evaluate the phase of the s-domain tranfer function for a given Complex value of s
1149 public double evalPhaseTransFunctS(Complex sValue){
1150 if(!this.padeAdded)this.transferPolesZeros();
1151 this.sValue = Complex.copy(sValue);
1152 Complex num = this.sNumer.evaluate(sValue);
1153 Complex den = this.sDenom.evaluate(sValue);
1154 Complex lagterm = Complex.plusOne();
1155 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1156 return (num.over(den).times(lagterm)).arg();
1157 }
1158
1159 // Evaluate the phase of the s-domain tranfer function for a sine wave input at a given frequency (s^-1)
1160 public double evalPhaseTransFunctS(double freq){
1161 if(!this.padeAdded)this.transferPolesZeros();
1162 this.sValue.reset(0.0D, 2.0D*Math.PI*freq);
1163 Complex num = this.sNumer.evaluate(this.sValue);
1164 Complex den = this.sDenom.evaluate(this.sValue);
1165 Complex lagterm = Complex.plusOne();
1166 if(this.deadTime!=0)lagterm = Complex.exp(this.sValue.times(-this.deadTime));
1167 return (num.over(den).times(lagterm)).arg();
1168 }
1169
1170 // Evaluate the z-domain tranfer function for the present value of z
1171 public Complex evalTransFunctZ(){
1172 Complex num = this.zNumer.evaluate(this.zValue);
1173 Complex den = this.zDenom.evaluate(this.zValue);
1174 return num.over(den);
1175 }
1176
1177 // Evaluate the z-domain tranfer function for a given Complex value of z
1178 public Complex evalTransFunctZ(Complex zValue){
1179 this.zValue = Complex.copy(zValue);
1180 Complex num = this.zNumer.evaluate(zValue);
1181 Complex den = this.zDenom.evaluate(zValue);
1182 return num.over(den);
1183 }
1184
1185 // Evaluate the magnitude of the z-domain tranfer function for the present value of z
1186 public double evalMagTransFunctZ(){
1187 Complex num = this.zNumer.evaluate(this.zValue);
1188 Complex den = this.zDenom.evaluate(this.zValue);
1189 return num.over(den).abs();
1190 }
1191
1192 // Evaluate the magnitude of the z-domain tranfer function for a given Complex value of z
1193 public double evalMagTransFunctZ(Complex zValue){
1194 this.zValue = Complex.copy(zValue);
1195 Complex num = this.zNumer.evaluate(zValue);
1196 Complex den = this.zDenom.evaluate(zValue);
1197 return num.over(den).abs();
1198 }
1199
1200 // Evaluate the phase of the z-domain tranfer function for the present value of z
1201 public double evalPhaseTransFunctZ(){
1202 Complex num = this.zNumer.evaluate(this.zValue);
1203 Complex den = this.zDenom.evaluate(this.zValue);
1204 return num.over(den).arg();
1205 }
1206
1207 // Evaluate the phase of the z-domain tranfer function for a given Complex value of z
1208 public double evalPhaseTransFunctZ(Complex zValue){
1209 this.zValue = Complex.copy(zValue);
1210 Complex num = this.zNumer.evaluate(zValue);
1211 Complex den = this.zDenom.evaluate(zValue);
1212 return num.over(den).arg();
1213 }
1214
1215 // Get the integration method option
1216 public int getIntegMethod(){
1217 return this.integMethod;
1218 }
1219
1220 // Get the z transform method option
1221 public int getZtransformMethod(){
1222 return this.ztransMethod;
1223 }
1224
1225 // Get the length of the time, input (time domain) and output (time domain) arrays
1226 public int getSampleLength(){
1227 return this.sampLen;
1228 }
1229
1230 // Get the forgetting factor
1231 public double getForgetFactor(){
1232 return this.forgetFactor;
1233 }
1234
1235 // Get the current time
1236 public double getCurrentTime(){
1237 return this.time[this.sampLen-1];
1238 }
1239
1240 // Get the time array
1241 public double[] getTime(){
1242 return this.time;
1243 }
1244
1245 // Get the current time domain input
1246 public double getCurrentInputT(){
1247 return this.inputT[this.sampLen-1];
1248 }
1249
1250 // Get the time domain input array
1251 public double[] getInputT(){
1252 return this.inputT;
1253 }
1254
1255 // Get the s-domain input
1256 public Complex getInputS(){
1257 return this.inputS;
1258 }
1259
1260 // Get the sampling period
1261 public double getDeltaT(){
1262 return this.deltaT;
1263 }
1264
1265 // Get the sampling frequency
1266 public double getSampFreq(){
1267 return this.sampFreq;
1268 }
1269
1270 // Get the Laplacian s value
1271 public Complex getS(){
1272 return this.sValue;
1273 }
1274
1275 // Get the z-transform z value
1276 public Complex getZ(){
1277 return this.zValue;
1278 }
1279
1280 // Get the degree of the original s-domain numerator polynomial
1281 public int getSnumerDeg(){
1282 return this.sNumerDeg;
1283 }
1284
1285 // Get the degree of the s-domain numerator polynomial after any dead time Pade approximation added
1286 public int getSnumerPadeDeg(){
1287 return this.sNumerDegPade;
1288 }
1289
1290 // Get the degree of the original s-domain denominator polynomial
1291 public int getSdenomDeg(){
1292 return this.sDenomDeg;
1293 }
1294
1295 // Get the degree of the s-domain denominator polynomial after any dead time Pade approximation added
1296 public int getSdenomPadeDeg(){
1297 return this.sDenomDegPade;
1298 }
1299
1300 // Get the original s-domain numerator polynomial
1301 public ComplexPoly getSnumer(){
1302 return this.sNumer.times(this.sNumerWorkingFactor);
1303 }
1304
1305 // Get the s-domain numerator polynomial after any dead time Pade approximation added
1306 public ComplexPoly getSnumerPade(){
1307 return this.sNumerPade.times(this.sNumerWorkingFactor);
1308 }
1309
1310 // Get the original s-domain denominator polynomial
1311 public ComplexPoly getSdenom(){
1312 return this.sDenom.times(this.sDenomWorkingFactor);
1313 }
1314
1315 // Get the s-domain denominator polynomial after any dead time Pade approximation added
1316 public ComplexPoly getSdenomPade(){
1317 return this.sDenomPade.times(this.sDenomWorkingFactor);
1318 }
1319
1320 // Get the degree of the z-domain numerator polynomial
1321 public int getZnumerDeg(){
1322 return this.zNumerDeg;
1323 }
1324
1325 // Get the degree of the z-domain denominator polynomial
1326 public int getZdenomDeg(){
1327 return this.zDenomDeg;
1328 }
1329
1330 // Get the z-domain numerator polynomial
1331 public ComplexPoly getZnumer(){
1332 return this.zNumer;
1333 }
1334
1335 // Get the z-domain denominator polynomial
1336 public ComplexPoly getZdenom(){
1337 return this.zDenom;
1338 }
1339
1340 // Get the s-domain zeros without any Pade zeros
1341 public Complex[] getZerosS(){
1342 if(this.sZeros==null)this.calcPolesZerosS();
1343 if(this.sZeros==null){
1344 System.out.println("Method BlackBox.getZerosS:");
1345 System.out.println("There are either no s-domain zeros for this transfer function");
1346 System.out.println("or the s-domain numerator polynomial has not been set");
1347 System.out.println("null returned");
1348 return null;
1349 }
1350 else{
1351 return this.sZeros;
1352 }
1353
1354 }
1355
1356 // Get the s-domain zeros plusany Pade zeros
1357 public Complex[] getZerosPadeS(){
1358 if(this.sZeros==null)this.calcPolesZerosS();
1359 if(!this.padeAdded)this.transferPolesZeros();
1360 if(this.sZerosPade==null){
1361 System.out.println("Method BlackBox.getZerosPadeS:");
1362 System.out.println("There are either no s-domain zeros for this transfer function");
1363 System.out.println("or the s-domain numerator polynomial has not been set");
1364 System.out.println("null returned");
1365 return null;
1366 }
1367 else{
1368 return this.sZerosPade;
1369 }
1370 }
1371
1372 // Get the s-domain poles without any Pade poles
1373 public Complex[] getPolesS(){
1374 if(this.sPoles==null)this.calcPolesZerosS();
1375 if(this.sPoles==null){
1376 System.out.println("Method BlackBox.getPolesS:");
1377 System.out.println("There are either no s-domain poles for this transfer function");
1378 System.out.println("or the s-domain denominator polynomial has not been set");
1379 System.out.println("null returned");
1380 return null;
1381 }
1382 else{
1383 return this.sPoles;
1384 }
1385 }
1386
1387 // Get the s-domain poles plus any Pade poles
1388 public Complex[] getPolesPadeS(){
1389 if(this.sPoles==null)this.calcPolesZerosS();
1390 if(!this.padeAdded)this.transferPolesZeros();
1391 if(this.sPolesPade==null){
1392 System.out.println("Method BlackBox.getPolesPadeS:");
1393 System.out.println("There are either no s-domain poles for this transfer function");
1394 System.out.println("or the s-domain denominator polynomial has not been set");
1395 System.out.println("null returned");
1396 return null;
1397 }
1398 else{
1399 return this.sPolesPade;
1400 }
1401 }
1402
1403
1404 // Get the z-domain zeros
1405 public Complex[] getZerosZ(){
1406 if(this.zZeros==null){
1407 System.out.println("Method BlackBox.getZerosZ:");
1408 System.out.println("There are either no z-domain zeros for this transfer function");
1409 System.out.println("or the z-domain numerator polynomial has not been set");
1410 System.out.println("null returned");
1411 return null;
1412 }
1413 else{
1414 return this.zZeros;
1415 }
1416 }
1417
1418 // Get the z-domain poles
1419 public Complex[] getPolesZ(){
1420 if(this.zPoles==null){
1421 System.out.println("Method BlackBox.getPolesZ:");
1422 System.out.println("There are either no z-domain poles for this transfer function");
1423 System.out.println("or the z-domain denominator polynomial has not been set");
1424 System.out.println("null returned");
1425 return null;
1426 }
1427 else{
1428 return this.zPoles;
1429 }
1430 }
1431
1432 // Get the map infinity zeros to zero or -1 option
1433 // maptozero: if true infinity s zeros map to zero
1434 // if false infinity s zeros map to minus one
1435 public boolean getMaptozero(){
1436 return this.maptozero;
1437 }
1438
1439 // Get the name of the black box
1440 public String getName(){
1441 return this.name;
1442 }
1443
1444 // Plot the poles and zeros of the BlackBox transfer function in the s-domain
1445 // Excludes any Pade poles and zeros
1446 public void plotPoleZeroS(){
1447 if(this.sNumer==null)throw new IllegalArgumentException("s domain numerator has not been set");
1448 if(this.sDenom==null)throw new IllegalArgumentException("s domain denominator has not been set");
1449 PlotPoleZero ppz = new PlotPoleZero(this.sNumer, this.sDenom);
1450 ppz.setS();
1451 ppz.pzPlot(this.name);
1452 }
1453
1454 // Plot the poles and zeros of the BlackBox transfer function in the s-domain
1455 // Includes Pade poles and zeros
1456 public void plotPoleZeroPadeS(){
1457 if(!this.padeAdded)this.transferPolesZeros();
1458 if(this.sNumerPade==null)throw new IllegalArgumentException("s domain numerator has not been set");
1459 if(this.sDenomPade==null)throw new IllegalArgumentException("s domain denominator has not been set");
1460 PlotPoleZero ppz = new PlotPoleZero(this.sNumerPade, this.sDenomPade);
1461 ppz.setS();
1462 ppz.pzPlot(this.name);
1463 }
1464
1465 // Plot the poles and zeros of the BlackBox transfer function in the z-domain
1466 public void plotPoleZeroZ(){
1467 PlotPoleZero ppz = new PlotPoleZero(this.zNumer, this.zDenom);
1468 if(this.zNumer==null)throw new IllegalArgumentException("z domain numerator has not been set");
1469 if(this.zDenom==null)throw new IllegalArgumentException("z domain denominator has not been set");
1470 ppz.setZ();
1471 ppz.pzPlot(this.name);
1472 }
1473
1474 // Bode plots for the magnitude and phase of the s-domain transfer function
1475 public void plotBode(double lowFreq, double highFreq){
1476 if(!this.padeAdded)this.transferPolesZeros();
1477 int nPoints = 100;
1478 double[][] cdata = new double[2][nPoints];
1479 double[] logFreqArray = new double[nPoints+1];
1480 double logLow = Fmath.log10(2.0D*Math.PI*lowFreq);
1481 double logHigh = Fmath.log10(2.0D*Math.PI*highFreq);
1482 double incr = (logHigh - logLow)/((double)nPoints-1.0D);
1483 double freqArray = lowFreq;
1484 logFreqArray[0]=logLow;
1485 for(int i=0; i<nPoints; i++){
1486 freqArray=Math.pow(10,logFreqArray[i]);
1487 cdata[0][i]=logFreqArray[i];
1488 cdata[1][i]=20.0D*Fmath.log10(this.evalMagTransFunctS(freqArray/(2.0*Math.PI)));
1489 logFreqArray[i+1]=logFreqArray[i]+incr;
1490 }
1491
1492 PlotGraph pgmag = new PlotGraph(cdata);
1493 pgmag.setGraphTitle("Bode Plot = magnitude versus log10[radial frequency]");
1494 pgmag.setGraphTitle2(this.name);
1495 pgmag.setXaxisLegend("Log10[radial frequency]");
1496 pgmag.setYaxisLegend("Magnitude[Transfer Function]");
1497 pgmag.setYaxisUnitsName("dB");
1498 pgmag.setPoint(0);
1499 pgmag.setLine(3);
1500 pgmag.plot();
1501 for(int i=0; i<nPoints; i++){
1502 freqArray=Math.pow(10,logFreqArray[i]);
1503 cdata[0][i]=logFreqArray[i];
1504 cdata[1][i]=this.evalPhaseTransFunctS(freqArray)*180.0D/Math.PI;
1505 }
1506 PlotGraph pgphase = new PlotGraph(cdata);
1507 pgphase.setGraphTitle("Bode Plot = phase versus log10[radial frequency]");
1508 pgphase.setGraphTitle2(this.name);
1509 pgphase.setXaxisLegend("Log10[radial frequency]");
1510 pgphase.setYaxisLegend("Phase[Transfer Function]");
1511 pgphase.setYaxisUnitsName("degrees");
1512 pgphase.setPoint(0);
1513 pgmag.setLine(3);
1514 pgphase.plot();
1515
1516 }
1517
1518 // Get the current time domain output for a given input and given time
1519 // resets deltaT
1520 public double getCurrentOutputT(double ttime, double inp){
1521 this.setInputT(ttime, inp);
1522 return this.getCurrentOutputT();
1523 }
1524
1525 // Get the current time domain output for the stored input
1526 public double getCurrentOutputT(){
1527 if(!this.padeAdded)this.transferPolesZeros();
1528
1529 ComplexPoly numerI = this.sNumerPade.times(new Complex(this.inputT[this.sampLen-1], 0.0));
1530 Complex[] polyC = {Complex.zero(), Complex.plusOne()};
1531 ComplexPoly polyH = new ComplexPoly(polyC);
1532 ComplexPoly denomI = this.sDenomPade.times(polyH);
1533
1534 Complex[][] coeffT = BlackBox.inverseTransform(numerI, denomI, this.sNumerWorkingFactor, this.sDenomScaleFactor);
1535
1536 Complex tempc = Complex.zero();
1537 for(int j=0; j<coeffT[0].length; j++){
1538 tempc.plusEquals(BlackBox.timeTerm(this.time[this.sampLen-1], coeffT[0][j], coeffT[1][j], coeffT[2][j]));
1539 }
1540 double outReal = tempc.getReal();
1541 double outImag = tempc.getImag();
1542 double temp;
1543 boolean outTest=true;
1544 if(outImag==0.0D)outTest=false;
1545 if(outTest){
1546 temp=Math.max(Math.abs(outReal),Math.abs(outImag));
1547 if(Math.abs((outReal-outImag)/temp)>1.e-5){
1548 outTest=false;
1549 }
1550 else{
1551 System.out.println("output in Blackbox.getCurrentOutputT() has a significant imaginary part");
1552 System.out.println("time = " + this.time[this.sampLen-1] + " real = " + outReal + " imag = " + outImag);
1553 System.out.println("Output equated to the real part");
1554 }
1555 }
1556 //for(int i=0; i<this.sampLen-2; i++)this.outputT[i]=this.outputT[i+1];
1557 this.outputT[this.sampLen-1] = outReal;
1558 return this.outputT[this.sampLen-1];
1559 }
1560
1561 // Get the time domain output array
1562 public double[] getOutputT(){
1563 return this.outputT;
1564 }
1565
1566 // Get the s-domain output for the stored input and s value.
1567 public Complex getOutputS(){
1568 if(!this.padeAdded)this.transferPolesZeros();
1569 Complex num = this.sNumer.evaluate(this.sValue);
1570 Complex den = this.sDenom.evaluate(this.sValue);
1571 this.outputS = num.over(den).times(this.inputS);
1572 if(this.deadTime!=0)this.outputS = this.outputS.times(Complex.exp(this.sValue.times(-this.deadTime)));
1573 return this.outputS;
1574 }
1575
1576 // Get the s-domain output for a given s value and input.
1577 public Complex getOutputS(Complex svalue, Complex inputs){
1578 if(!this.padeAdded)this.transferPolesZeros();
1579 this.inputS = inputs;
1580 this.sValue = svalue;
1581 Complex num = this.sNumer.evaluate(this.sValue);
1582 Complex den = this.sDenom.evaluate(this.sValue);
1583 this.outputS = num.over(den).times(this.inputS);
1584 if(this.deadTime!=0)this.outputS = this.outputS.times(Complex.exp(this.sValue.times(-this.deadTime)));
1585 return this.outputS;
1586 }
1587
1588 // Reset the number of points used in plotting a response curve
1589 public void setNplotPoints(int nPoints){
1590 this.nPlotPoints = nPoints;
1591 }
1592
1593 // Return the number of points used in plotting a response curve
1594 public int getNplotPoints(){
1595 return this.nPlotPoints;
1596 }
1597
1598 // Plots the time course for an impulse input
1599 public void impulseInput(double impulseMag, double finalTime){
1600 if(!this.padeAdded)this.transferPolesZeros();
1601
1602 // Multiply transfer function by impulse magnitude (impulseMag)
1603 ComplexPoly impulseN = new ComplexPoly(0);
1604 impulseN.resetCoeff(0, Complex.plusOne().times(impulseMag));
1605 ComplexPoly numerT = this.sNumerPade.times(impulseN);
1606 ComplexPoly denomT = this.sDenomPade.copy();
1607 String graphtitle1 = "Impulse Input Transient: Impulse magnitude = "+impulseMag;
1608 String graphtitle2 = this.getName();
1609 BlackBox.transientResponse(this.nPlotPoints, finalTime, this.deadTime, numerT, denomT, graphtitle1, graphtitle2, this.sNumerWorkingFactor, this.sDenomScaleFactor);
1610 }
1611
1612 // Plots the time course for a unit impulse input
1613 public void impulseInput(double finalTime){
1614 this.impulseInput(1.0D, finalTime);
1615 }
1616
1617 // Plots the time course for a step input
1618 public void stepInput(double stepMag, double finalTime){
1619 Complex sNumer0 = this.sNumerPade.coeffCopy(0);
1620 Complex sDenom0 = this.sDenomPade.coeffCopy(0);
1621 boolean test0 = false;
1622 if(Complex.isReal(sNumer0) && Complex.isReal(sDenom0))test0=true;
1623
1624 if(sNumerDeg==0 && sDenomDeg==0 && test0){
1625 // Calculate time course outputs
1626 int n = 51; // number of points on plot
1627 double incrT = finalTime/(double)(n-2); // plotting increment
1628 double cdata[][] = new double [2][n]; // plotting array
1629
1630 cdata[0][0]=0.0D;
1631 cdata[0][1]=0.0D;
1632 for(int i=2; i<n; i++){
1633 cdata[0][i]=cdata[0][i-1]+incrT;
1634 }
1635 double kpterm = sNumer0.getReal()*stepMag/sDenom0.getReal();
1636 cdata[1][0]=0.0D;
1637 for(int i=1; i<n; i++){
1638 cdata[1][i] = kpterm;
1639 }
1640 if(this.deadTime!=0.0D)for(int i=0; i<n; i++)cdata[0][i] += this.deadTime;
1641
1642 // Plot
1643 PlotGraph pg = new PlotGraph(cdata);
1644
1645 pg.setGraphTitle("Step Input Transient: Step magnitude = "+stepMag);
1646 pg.setGraphTitle2(this.getName());
1647 pg.setXaxisLegend("Time");
1648 pg.setXaxisUnitsName("s");
1649 pg.setYaxisLegend("Output");
1650 pg.setPoint(0);
1651 pg.setLine(3);
1652 pg.plot();
1653
1654 }
1655 else{
1656 if(!this.padeAdded)this.transferPolesZeros();
1657 // Multiply transfer function by step magnitude (stepMag)/s
1658 ComplexPoly numerT = this.sNumer.times(stepMag);
1659 Complex[] polyC = {Complex.zero(), Complex.plusOne()};
1660 ComplexPoly polyH = new ComplexPoly(polyC);
1661 ComplexPoly denomT = this.sDenom.times(polyH);
1662 String graphtitle1 = "Step Input Transient: Step magnitude = "+stepMag;
1663 String graphtitle2 = this.getName();
1664
1665 BlackBox.transientResponse(this.nPlotPoints, finalTime, this.deadTime, numerT, denomT, graphtitle1, graphtitle2, this.sNumerWorkingFactor, this.sDenomScaleFactor);
1666 }
1667 }
1668
1669 // Plots the time course for a unit step input
1670 public void stepInput(double finalTime){
1671 this.stepInput(1.0D, finalTime);
1672 }
1673
1674 // Plots the time course for an nth order ramp input (a.t^n)
1675 public void rampInput(double rampGradient, int rampOrder, double finalTime){
1676 if(!this.padeAdded)this.transferPolesZeros();
1677
1678 // Multiply transfer function by ramp input (rampGradient)(rampOrder!)/s^(ramporder+1)
1679 ComplexPoly numerT = this.sNumer.times(rampGradient*Fmath.factorial(rampOrder));
1680 Complex[] polyC = Complex.oneDarray(rampOrder+1);
1681 for(int i=0; i<rampOrder; i++)polyC[i] = Complex.zero();
1682 polyC[rampOrder] = Complex.plusOne();
1683 ComplexPoly polyH = new ComplexPoly(polyC);
1684 ComplexPoly denomT = this.sDenom.times(polyH);
1685 String graphtitle1 = "";
1686 if(rampGradient!=1.0D){
1687 if(rampOrder!=1){
1688 graphtitle1 += "nth order ramp (at^n) input transient: a = "+rampGradient+" n = "+rampOrder;
1689 }
1690 else{
1691 graphtitle1 += "First order ramp (at) input transient: a = "+rampGradient;
1692 }
1693 }
1694 else{
1695 if(rampOrder!=1){
1696 graphtitle1 += "Unit ramp (t) input transient";
1697 }
1698 else{
1699 graphtitle1 += "nth order ramp (t^n) input transient: n = "+rampOrder;
1700 }
1701 }
1702 String graphtitle2 = this.getName();
1703 BlackBox.transientResponse(this.nPlotPoints, finalTime, this.deadTime, numerT, denomT, graphtitle1, graphtitle2, this.sNumerWorkingFactor, this.sDenomScaleFactor);
1704 }
1705
1706 // Plots the time course for an nth order ramp input (t^n)
1707 public void rampInput(int rampOrder, double finalTime){
1708 double rampGradient = 1.0D;
1709 this.rampInput(rampGradient, rampOrder, finalTime);
1710 }
1711
1712 // Plots the time course for a first order ramp input (at)
1713 public void rampInput(double rampGradient, double finalTime){
1714 int rampOrder = 1;
1715 this.rampInput(rampGradient, rampOrder, finalTime);
1716 }
1717
1718 // Plots the time course for a unit ramp input (t)
1719 public void rampInput(double finalTime){
1720 double rampGradient = 1.0D;
1721 int rampOrder = 1;
1722 this.rampInput(rampGradient, rampOrder, finalTime);
1723 }
1724
1725 // Plots the time course for a given transfer function from time t = zero for a quiescent system
1726 // Denominator scaling factor calculated
1727 public static void transientResponse(int nPoints, double finalTime, double deadTime, ComplexPoly numerT, ComplexPoly denomT, String graphtitle1, String graphtitle2){
1728 Complex[] roots = denomT.rootsNoMessages();
1729 Complex magDenom = BlackBox.scaleFactor(denomT, roots);
1730 Complex magNumer = Complex.plusOne();
1731 BlackBox.transientResponse(nPoints, finalTime, deadTime, numerT, denomT, graphtitle1, graphtitle2, magNumer, magDenom);
1732 }
1733
1734 // Plots the time course for a given transfer function from time t = zero for a quiescent system
1735 // Denominator scaling factor provided
1736 public static void transientResponse(int nPoints, double finalTime, double deadTime, ComplexPoly numerT, ComplexPoly denomT, String graphtitle1, String graphtitle2, Complex magN, Complex magD){
1737 // Obtain coefficients and constants of an partial fraction expansion
1738
1739
1740 Complex[][] coeffT = BlackBox.inverseTransform(numerT, denomT, magN, magD);
1741
1742 // Calculate time course outputs
1743 int m = denomT.getDeg(); // number of Aexp(-at) terms
1744 double incrT = finalTime/(double)(nPoints-1); // plotting increment
1745 double cdata[][] = new double [2][nPoints]; // plotting array
1746 double temp = 0.0D; // working variable
1747 Complex tempc = new Complex(); // working variable
1748 double outReal = 0.0D; // real part of output
1749 double outImag = 0.0D; // imaginary part of output (should be zero)
1750 boolean outTest = true; // false if outImag=zero
1751
1752 cdata[0][0]=0.0D;
1753 for(int i=1; i<nPoints; i++){
1754 cdata[0][i]=cdata[0][i-1]+incrT;
1755 }
1756 for(int i=0; i<nPoints; i++){
1757 outTest= true;
1758 tempc = Complex.zero();
1759 for(int j=0; j<m; j++){
1760 tempc.plusEquals(BlackBox.timeTerm(cdata[0][i], coeffT[0][j], coeffT[1][j], coeffT[2][j]));
1761 }
1762 outReal = tempc.getReal();
1763 outImag = tempc.getImag();
1764 if(outImag==0.0D)outTest=false;
1765 if(outTest){
1766 temp=Math.max(Math.abs(outReal),Math.abs(outImag));
1767 if(Math.abs((outReal-outImag)/temp)>1.e-5){
1768 outTest=false;
1769 }
1770 else{
1771 System.out.println("output in Blackbox.stepInput has a significant imaginary part");
1772 System.out.println("time = " + cdata[0][i] + " real = " + outReal + " imag = " + outImag);
1773 System.out.println("Output equated to the real part");
1774 }
1775 }
1776 cdata[1][i]=outReal;
1777 cdata[0][i]+=deadTime;
1778 }
1779
1780 // Plot
1781 PlotGraph pg = new PlotGraph(cdata);
1782
1783 pg.setGraphTitle(graphtitle1);
1784 pg.setGraphTitle2(graphtitle2);
1785 pg.setXaxisLegend("Time");
1786 pg.setXaxisUnitsName("s");
1787 pg.setYaxisLegend("Output");
1788 pg.setPoint(0);
1789 pg.setLine(3);
1790 pg.setNoYoffset(true);
1791 if(deadTime<(cdata[0][nPoints-1]-cdata[0][0]))pg.setNoXoffset(true);
1792 pg.setXlowFac(0.0D);
1793 pg.setYlowFac(0.0D);
1794 pg.plot();
1795 }
1796
1797
1798 // Returns the output term for a given time, coefficient, constant and power
1799 // for output = A.time^(n-1).exp(constant*time)/(n-1)!
1800 // Complex arguments and return
1801 public static Complex timeTerm(double ttime, Complex coeff, Complex constant, Complex power){
1802 Complex ret = new Complex();
1803 int n = (int)power.getReal() - 1;
1804 ret = coeff.times(Math.pow(ttime,n));
1805 ret = ret.over(Fmath.factorial(n));
1806 ret = ret.times(Complex.exp(constant.times(ttime)));
1807 return ret;
1808 }
1809
1810 // Returns the output term for a given time, coefficient, constant and power
1811 // for output = A.time^(n-1).exp(constant*time)/(n-1)!
1812 // Real arguments and return
1813 public static double timeTerm(double ttime, double coeff, double constant, int power){
1814 int n = power - 1;
1815 double ret = coeff*Math.pow(ttime,n);
1816 ret = ret/Fmath.factorial(n);
1817 ret = ret*(Math.exp(constant*ttime));
1818 return ret;
1819 }
1820
1821 // Returns the output term for a given time, coefficient, constant and power
1822 // for output = A.time^(n-1).exp(constant*time)/(n-1)!
1823 // Real arguments and return - all double
1824 public static double timeTerm(double ttime, double coeff, double constant, double power){
1825 double n = power - 1;
1826 double ret = coeff*Math.pow(ttime,n);
1827 ret = ret/Fmath.factorial(n);
1828 ret = ret*(Math.exp(constant*ttime));
1829 return ret;
1830 }
1831
1832 // Returns the coefficients A, the constant a and the power n in the f(A.exp(-at),n) term for the
1833 // the inverse Laplace transform of a complex polynolial divided
1834 // by a complex polynomial expanded as partial fractions
1835 // A and a are returnd as a 2 x n double array were n is the number of terms
1836 // in the partial fraction. the first row contains the A values, the second the a values
1837 // denominator scaling factor calculated
1838 public static double[][] inverseTransformToReal(ComplexPoly numer, ComplexPoly denom){
1839 Complex[][] com = inverseTransform(numer, denom);
1840 int n = com[0].length;
1841 double[][] ret = new double[3][n];
1842 for(int i=0; i<n;i++){
1843 ret[0][i] = com[0][i].getReal();
1844 if(Math.abs((ret[0][i]-com[0][i].getImag())/ret[0][i])>1.e-5){
1845 System.out.println("BlackBox inverseTransformToReal coefficient A[" + i + "] has a significant imaginary part: " + com[0][i]);
1846 System.out.println("A equated to the real part");
1847 System.out.println("inverseTransform method may be more appropriate");
1848 }
1849 ret[1][i] = com[1][i].getReal();
1850 if(Math.abs((ret[1][i]-com[1][i].getImag())/ret[1][i])>1.e-5){
1851 System.out.println("BlackBox inverseTransformToReal coefficient a[" + i + "] has a significant imaginary part: " + com[1][i]);
1852 System.out.println("a equated to the real part");
1853 System.out.println("inverseTransform method may be more appropriate");
1854 }
1855 ret[2][i] = com[2][i].getReal();
1856 }
1857 return ret;
1858 }
1859
1860
1861
1862 // Returns the coefficients A, the constant a and the power n in the f(A.exp(-at),n) term for the
1863 // the inverse Laplace transform of a complex polynolial divided
1864 // by a complex polynomial expanded as partial fractions
1865 // A and a are returnd as a 3 x n Complex array were n is the number of terms
1866 // in the partial fraction. The first row contains the A values, the second the a values and the third the power of the denominator
1867 // denominator scaling factor calculated
1868 public static Complex[][] inverseTransform(ComplexPoly numer, ComplexPoly denom){
1869 Complex[] roots = denom.rootsNoMessages();
1870 Complex magDenom = BlackBox.scaleFactor(denom, roots);
1871 Complex magNumer = Complex.plusOne();
1872 return inverseTransform(numer, denom, magNumer, magDenom);
1873 }
1874
1875 // Returns the coefficients A, the constant a and the power n in the f(A.exp(-at),n) term for the
1876 // the inverse Laplace transform of a complex polynolial divided
1877 // by a complex polynomial expanded as partial fractions
1878 // A and a are returnd as a 3 x n Complex array were n is the number of terms
1879 // in the partial fraction. The first row contains the A values, the second the a values and the third the power of the denominator
1880 // denominator scaling factor provided
1881 public static Complex[][] inverseTransform(ComplexPoly numer, ComplexPoly denom, Complex magNumer, Complex magDenom){
1882
1883 int polesN = denom.getDeg(); // number of poles
1884 int zerosN = numer.getDeg(); // numer of zeros
1885 if(zerosN>=polesN)throw new IllegalArgumentException("The degree of the numerator is equal to or greater than the degree of the denominator");
1886 Complex[][] ret = Complex.twoDarray(3, polesN); // array for returning coefficients, constants and powers
1887
1888 // Special case: input = A/(B + C)s
1889 if(polesN==1 && zerosN==0){
1890 Complex num = numer.coeffCopy(0);
1891 Complex den0 = denom.coeffCopy(0);
1892 Complex den1 = denom.coeffCopy(1);
1893 ret[0][0] = num.over(den1);
1894 ret[1][0] = Complex.minusOne().times(den0.over(den1));
1895 ret[2][0] = new Complex(1.0, 0.0);
1896 return ret;
1897 }
1898
1899 int nDifferentRoots = polesN; // number of roots of different values
1900 int nSetsIdenticalRoots = 0; // number of sets roots of identical value
1901 Complex[] poles = denom.rootsNoMessages(); // poles array
1902 int[] polePower = new int[polesN]; // power, n, of each (s - root)^n term
1903 boolean[] poleSet = new boolean[polesN]; // true if root has been identified as either equal to another root
1904 int[] poleIdent = new int[polesN]; // same integer for identical (s-root) terms; integer = index of first case of that root
1905 int[] poleHighestPower = new int[polesN]; // highest pole power for that set of identical poles
1906 boolean[] termSet = new boolean[polesN]; // false if n in (s-root)^n is greater than 1 and less than maximum value of n for that root
1907 double identicalRootLimit = 1.0e-2; // roots treated as identical if equal to one part in identicalRootLimit
1908 int[] numberInSet = new int[polesN]; // number of poles indentical to this pole including this pole
1909
1910 // Find identical roots within identicalRootLimit and assign power n [ (s-a)^n] to all roots
1911 int power = 0;
1912 Complex identPoleAverage = new Complex();
1913 int lastPowerIndex=0;
1914 for(int i=0; i<polesN; i++)poleSet[i]=false;
1915 for(int i=0; i<polesN; i++)termSet[i]=true;
1916 for(int i=0; i<polesN; i++){
1917 if(!poleSet[i]){
1918 power=1;
1919 polePower[i]=1;
1920 poleHighestPower[i]= 1;
1921 poleIdent[i]=i;
1922 numberInSet[i]=1;
1923 identPoleAverage = poles[i];
1924 for(int j=i+1; j<polesN; j++){
1925 if(!poleSet[j]){
1926 if(poles[i].isEqualWithinLimits(poles[j],identicalRootLimit)){
1927 poleIdent[j]=i;
1928 polePower[j]=++power;
1929 poleSet[j]=true;
1930 poleSet[i]=true;
1931 termSet[j]=false;
1932 termSet[i]=false;
1933 lastPowerIndex=j;
1934 nDifferentRoots--;
1935 identPoleAverage = identPoleAverage.plus(poles[j]);
1936 }
1937 else{
1938 poleIdent[j]=j;
1939 polePower[j]=1;
1940 }
1941 }
1942 }
1943 }
1944
1945 if(poleSet[i]){
1946 nDifferentRoots--;
1947 nSetsIdenticalRoots++;
1948
1949 // Set termSet to true if pole is recurring term with the highest power
1950 termSet[lastPowerIndex]=true;
1951
1952 // Replace roots within identicalRootLimit with their average value
1953 identPoleAverage = identPoleAverage.over(power);
1954 for(int j=0; j<polesN; j++){
1955 if(poleSet[j] && poleIdent[j]==i){
1956 poles[j] = identPoleAverage;
1957 poleHighestPower[i] = power;
1958 numberInSet[j] = power;
1959 }
1960 }
1961 }
1962 }
1963
1964 // Calculate pole average
1965 Complex poleAverage = Complex.zero();
1966 Complex absPoleAverage = Complex.zero();
1967 for(int i=0; i<polesN; i++){
1968 poleAverage = poleAverage.plus(poles[i]);
1969 absPoleAverage = absPoleAverage.plus(poles[i].abs());
1970 }
1971 poleAverage = poleAverage.over(polesN);
1972 absPoleAverage = absPoleAverage.over(polesN);
1973
1974 // Calculate pole substitute for identical substitution values
1975 Complex poleSubstitute = poleAverage;
1976 if(poleSubstitute.isZero())poleSubstitute = absPoleAverage;
1977 if(poleSubstitute.isZero())poleSubstitute = Complex.plusOne();
1978
1979 // Choose initial set of s substitution values
1980 Complex[] subValues = Complex.oneDarray(polesN);
1981 boolean[] subSet = new boolean[polesN];
1982 for(int i=0; i<polesN; i++)subSet[i] = false;
1983
1984 Complex[] shifts = null;
1985 Complex delta = new Complex(1.7, 0.0); // root separation factor
1986 for(int i=0; i<polesN; i++)subValues[i] = poles[i].copy();
1987 int currentNumberInSet = 0;
1988 if(nSetsIdenticalRoots>0){
1989 for(int i=0; i<polesN; i++){
1990 if(numberInSet[i]>1 && !subSet[i]){
1991 currentNumberInSet = numberInSet[i];
1992 shifts = Complex.oneDarray(numberInSet[i]);
1993 int centre = numberInSet[i]/2;
1994 if(Fmath.isEven(numberInSet[i])){
1995 for(int j=0; j<centre; j++){
1996 shifts[centre+j] = delta.times((double)(j+1));
1997 shifts[centre-1-j] = shifts[centre+j].times(-1.0);
1998 }
1999 }
2000 else{
2001 shifts[centre] = Complex.zero();
2002 for(int j=0; j<centre; j++){
2003 shifts[centre+1+j] = delta.times((double)(j+1));
2004 shifts[centre-1-j] = shifts[centre+j].times(-1.0);
2005 }
2006 }
2007 int kk = 0;
2008 for(int j=0; j<polesN; j++){
2009 if(!subSet[j] && numberInSet[j]==currentNumberInSet){
2010 Complex incr = poles[j];
2011 if(incr.isZero())incr = poleSubstitute;
2012 subValues[j] = shifts[kk].times(incr);
2013 subSet[j] = true;
2014 kk++;
2015
2016 }
2017 }
2018 }
2019 }
2020 }
2021
2022 // Check for identical and very close substitution values
2023 boolean testii = true;
2024 int ii = 0;
2025 int nAttempts = 0;
2026 while(testii){
2027 int jj = ii + 1;
2028 boolean testjj = true;
2029 while(testjj){
2030 if(subValues[ii].isEqualWithinLimits(subValues[jj],identicalRootLimit)){
2031 subValues[ii] = subValues[ii].plus(poleSubstitute.times((double)nAttempts));
2032 nAttempts++;
2033 ii=0;
2034 testjj = false;
2035 if(nAttempts>1000000)throw new IllegalArgumentException("a non repeating set of substitution values could not be foumd");
2036 }
2037 else{
2038 jj++;
2039 }
2040 if(jj>=polesN)testjj = false;
2041 }
2042 ii++;
2043 if(ii>=polesN-1)testii = false;
2044 }
2045
2046 // Set up the linear equations
2047 // Create vector and matrix arrays
2048 Complex[][] mat = Complex.twoDarray(polesN, polesN);
2049 Complex[] vec = Complex.oneDarray(polesN);
2050
2051 // Fill vector
2052 for(int i=0; i<polesN; i++){
2053 if(zerosN>0){
2054 vec[i] = numer.evaluate(subValues[i]);
2055 }
2056 else{
2057 vec[i] = numer.coeffCopy(0);
2058 }
2059 }
2060
2061 // fill matrix
2062 for(int i=0; i<polesN; i++){
2063 for(int j=0; j<polesN; j++){
2064 Complex denomTerm = Complex.plusOne();
2065 int powerD = 0;
2066 for(int k=0; k<polesN; k++){
2067 if(termSet[k]){
2068 if(j!=k){
2069 if(polePower[k]==1){
2070 denomTerm = denomTerm.times(subValues[i].minus(poles[k]));
2071 }
2072 else{
2073 denomTerm = denomTerm.times(Complex.pow(subValues[i].minus(poles[k]), polePower[k]));
2074 }
2075 }
2076 else{
2077 if(polePower[j]<poleHighestPower[j]){
2078 powerD = poleHighestPower[j] - polePower[j];
2079 if(powerD==1){
2080 denomTerm = denomTerm.times(subValues[i].minus(poles[k]));
2081 }
2082 else{
2083 if(powerD!=0){
2084 denomTerm = denomTerm.times(Complex.pow(subValues[i].minus(poles[k]), powerD));
2085 }
2086 }
2087 }
2088 }
2089 }
2090 }
2091 mat[i][j] = denomTerm;
2092 }
2093
2094 }
2095
2096
2097 // Solve linear equations
2098 ComplexMatrix cmat = new ComplexMatrix(mat);
2099 Complex[] terms = cmat.solveLinearSet(vec);
2100
2101 // fill ret for returning
2102 for(int i=0; i<polesN; i++){
2103 ret[0][i]=terms[i].times(magNumer).over(magDenom);
2104 ret[1][i]=poles[i];
2105 ret[2][i].reset(polePower[i],0.0D);
2106 }
2107 return ret;
2108
2109 }
2110
2111 // Deep copy
2112 public BlackBox copy(){
2113 if(this==null){
2114 return null;
2115 }
2116 else{
2117 BlackBox bb = new BlackBox();
2118 this.copyBBvariables(bb);
2119 return bb;
2120 }
2121 }
2122
2123 // Copies BlackBox variables
2124 public void copyBBvariables(BlackBox bb){
2125
2126 bb.sampLen = this.sampLen;
2127 bb.inputT = Conv.copy(this.inputT);
2128 bb.outputT = Conv.copy(this.outputT);
2129 bb.time = Conv.copy(this.time);
2130 bb.forgetFactor = this.forgetFactor;
2131 bb.deltaT = this.deltaT;
2132 bb.sampFreq = this.sampFreq;
2133 bb.inputS = this.inputS.copy();
2134 bb.outputS = this.outputS.copy();
2135 bb.sValue = this.sValue.copy();
2136 bb.zValue = this.zValue.copy();
2137 bb.sNumer = this.sNumer.copy();
2138 bb.sDenom = this.sDenom.copy();
2139 bb.zNumer = this.zNumer.copy();
2140 bb.zDenom = this.zDenom.copy();
2141 bb.sNumerSet = this.sNumerSet;
2142 bb.sDenomSet = this.sDenomSet;
2143 bb.sNumerScaleFactor = this.sNumerScaleFactor;
2144 bb.sDenomScaleFactor = this.sDenomScaleFactor;
2145 bb.sPoles = Complex.copy(this.sPoles);
2146 bb.sZeros = Complex.copy(this.sZeros);
2147 bb.zPoles = Complex.copy(this.zPoles);
2148 bb.zZeros = Complex.copy(this.zZeros);
2149 bb.sNumerDeg = this.sNumerDeg;
2150 bb.sDenomDeg = this.sDenomDeg;
2151 bb.zNumerDeg = this.zNumerDeg;
2152 bb.zDenomDeg = this.zDenomDeg;
2153 bb.deadTime = this.deadTime;
2154 bb.orderPade = this.orderPade;
2155 bb.sNumerPade = this.sNumerPade.copy();
2156 bb.sDenomPade = this.sDenomPade.copy();
2157 bb.sPolesPade = Complex.copy(this.sPolesPade);
2158 bb.sZerosPade = Complex.copy(this.sZerosPade);
2159 bb.sNumerDegPade = this.sNumerDegPade;
2160 bb.sDenomDegPade = this.sDenomDegPade;
2161 bb.maptozero = this.maptozero;
2162 bb.padeAdded = this.padeAdded;
2163 bb.integrationSum = this.integrationSum;
2164 bb.integMethod = this.integMethod;
2165 bb.ztransMethod = this.ztransMethod;
2166 bb.name = this.name;
2167 bb.fixedName = this.fixedName;
2168 bb.nPlotPoints = this.nPlotPoints;
2169
2170 }
2171
2172
2173 // Clone - overrides Java.Object method clone
2174 public Object clone(){
2175 return (Object)this.copy();
2176 }
2177}
2178
Note: See TracBrowser for help on using the repository browser.