se.lth.control.algorithms

Class Polynomial

java.lang.Object

se.lth.control.algorithms.Polynomial

public class Polynomial
extends java.lang.Object

A simple polynomial package designed for use in control applications. The internal representation of a polynomial is a vector where the coefficient of the highest power is stored at the first position, i.e. position zero. A polynomial

  A=a_0 z^n + a_1 z^(n-1) + ... + a_n

  will thus be represented as

  A=[a_0, a_1, ... ,a_n]
  

The data type of the polynomial coefficients is double.

Constructor Summary

Constructors

Constructor and Description
Polynomial() Constructs a polynomial of degree zero with the only coefficient equal to zero.
Polynomial(int Degree, double[] coefficients) Constructs a polynomial of degree Degree with the coefficients from Coefficients.
Polynomial(int Degree, double[] coefficients, java.lang.String str) This constructor works in the same way as the one above, but has in addition the ability to give a name to the polynomial.
Polynomial(Polynomial p) Constructs a polynomial identical to the argument p.

Method Summary

Modifier and Type	Method and Description
Polynomial	add(Polynomial p2) Returns a Polynomial whose value is the quotient of (this + p2).
static Polynomial	add(Polynomial p1, Polynomial p2) Returns a Polynomial whose value is the quotient of (p1 + p2).
Polynomial	assign(Polynomial p) Gives the polynial the value of p.
java.lang.Object	clone() Creates an new identical Polynomial object.
void	cut(double reps) Normalizes the polynomial, i.e.
static int	diophantineMDS(Polynomial a, Polynomial b, Polynomial c, Polynomial xx, Polynomial yy) DiophantineMDS finds the minimal degree solution to the equation AR+BS=Ac.
void	display() Displays the polynomial on the console.
void	display(java.lang.String str) Displays the polynomial on the console.
Polynomial	div(double r) Returns a Polynomial whose value is the quotient of (this / r).
Polynomial	div(Polynomial p2) Returns a Polynomial whose value is the quotient of (this / p2).
static Polynomial	div(Polynomial p, double r) Returns a Polynomial whose value is (p / r).
static Polynomial	div(Polynomial p1, Polynomial p2) Returns a Polynomial whose value is the quotient of (p1 / p2).
static void	dyadicReduction(double[] a, double[] b, double[] alpha, double[] beta, int i0, int i1, int i2) Given vectors
boolean	equals(Polynomial p2) Compares this Polynomial with the specified Polynomial for equality.
double	eval(double r) Calculates the value of the polynomial for q = r, i.e.
static void	GCD(Polynomial a, Polynomial b, double reps, Polynomial g, Polynomial x, Polynomial y, Polynomial u, Polynomial v) This function calculates the greatest common divisor of polynomials A and B using Euclid's algorithm.
double	get(int elem) Returns the value of element i.
int	getDegree() Returns the degree of the polynomial.
static void	integralDesign(Polynomial A, Polynomial B, Polynomial Ac, Polynomial R, Polynomial S, double x0)
boolean	isStable() This method uses Jury's Criterion to determine if the polynomial is stable.
static void	LDFilter(double[] l, double[] d, double[] phi, double[] theta, double lambda, int degree) LFilter s an implementation of an recursive estimator and would together with any of the design functions above form a complete adaptive controller.
static void	LQGDesign(Polynomial A, Polynomial B, Polynomial C, Polynomial R, Polynomial S, double rho) Calculates a LQG-controller for the system $A$, $B$, and $C$ with the loss function coefficient $\rho$.
Polynomial	makeMonic() Makes the polynomial monic by dividing all coefficients with the highest power coefficient.
static int	MDPPNZCDesign(Polynomial A, Polynomial B, Polynomial Am, Polynomial Ao, Polynomial R, Polynomial S, Polynomial T) MDPPNZCDesign stands for Minimal Degree Pole Placement with No Zero Cancelation.
Polynomial	mod(Polynomial p2) Returns a Polynomial whose value is the remainder of (this / p2).
static Polynomial	mod(Polynomial p1, Polynomial p2) Returns a Polynomial whose value is the remainder of (p1 / p2).
Polynomial	mul(double r) Returns a Polynomial whose value is (this * r).
static Polynomial	mul(double r, Polynomial p) Returns a Polynomial whose value is (r * p).
Polynomial	mul(Polynomial p2) Returns a Polynomial whose value is (this * p2).
static Polynomial	mul(Polynomial p, double r) Returns a Polynomial whose value is (p * r).
static Polynomial	mul(Polynomial p1, Polynomial p2) Returns a Polynomial whose value is (p1 * p2).
Polynomial	reciprocal() Returns the reciprocal polynomial.
static void	robustDesign(Polynomial A, Polynomial B, Polynomial Ac, Polynomial R, Polynomial S, double x0)
static void	robustIntegralDesign(Polynomial A, Polynomial B, Polynomial Ac, Polynomial R, Polynomial S, double x1, double x0)
void	set(int elem, double value) Sets the value of element i.
void	set(int newDeg, double[] newCoeff) Sets the degree of the Polynomial to newDeg and uses the values in newCoeff as coefficients.
void	set(int Degree, double[] vals, java.lang.String str) Set the degree of the Polyomial to Degree, the coefficients vals, and the name to str.
void	setName(java.lang.String str) Sets the class attribute name to str.
static void	sfactorize(Polynomial b, Polynomial a) The reciprocal polynomial of A, denoted A* is obtained by reversing the order of the coefficients of A.
void	shiftBackward(double newFirstcoeff) This method shifts the coefficients of the polynomial to the right and sets the highest power coefficient to r.\\
void	shiftForward(double newLastcoeff) This method shifts the coefficients of the polynomial to the left and sets the last coefficient to r.
Polynomial	sub(Polynomial p2) Returns a Polynomial whose value is the quotient of (this - p2).
static Polynomial	sub(Polynomial p1, Polynomial p2) Returns a Polynomial whose value is the quotient of (p1 + p2).
java.lang.String	toString() Returns a string representation of the Polynomial

Methods inherited from class java.lang.Object

equals, finalize, getClass, hashCode, notify, notifyAll, wait, wait, wait

Constructor Detail

Polynomial

public Polynomial()

Constructs a polynomial of degree zero with the only coefficient equal to zero.

Polynomial

public Polynomial(Polynomial p)

Constructs a polynomial identical to the argument p.

Polynomial

public Polynomial(int Degree,
                  double[] coefficients)

Constructs a polynomial of degree Degree with the coefficients from Coefficients. The first element in both the input vector Coefficients and the resulting polynomial represents the highest power coefficient.

Polynomial

public Polynomial(int Degree,
                  double[] coefficients,
                  java.lang.String str)

This constructor works in the same way as the one above, but has in addition the ability to give a name to the polynomial. This feature is mainly used for debug purposes. The program example below uses this constructor to create the polynomial A=q^2+1.5q-2.

     
     int   deg;
     double[] coeff = new double[3];
     Polynomial A;
     
     deg = 2;
     coeff[0] = 1;
     coeff[1] = 1.5;
     coeff[2] = -2;
     
     A = new Polynomial(deg, coeff, "A");
     A.display();
         

     This program generates the following output.
     A = 1.0 q^2 + 1.5 q - 2.0
     

Method Detail

clone

public java.lang.Object clone()

Creates an new identical Polynomial object.

Overrides:

clone in class java.lang.Object

set

public void set(int Degree,
                double[] vals,
                java.lang.String str)

Set the degree of the Polyomial to Degree, the coefficients vals, and the name to str.

setName

public void setName(java.lang.String str)

Sets the class attribute name to str.

display

public void display()

Displays the polynomial on the console. For example, let the variable A be defined and make the following function call:

     A.setName("A");
     A.display();
     
     The following output is generated:
     
     A = 1.0 q^2 + 1.5 q - 2.0
     

display

public void display(java.lang.String str)

Displays the polynomial on the console. If the function is called with an argument it will write the string str followed by the polynomial to the console. For example, let the variable A be defined as in the program example above, and make the following function call:

     A.display("This is polynomial A:");
     
     The following output is generated:
     
     This is polynomial A:1.0 q^2 + 1.5 q - 2.0
     

toString

public java.lang.String toString()

Returns a string representation of the Polynomial

Overrides:

toString in class java.lang.Object

reciprocal

public Polynomial reciprocal()

Returns the reciprocal polynomial. The reciprocal polynomial is defined as follows:

     Let A=3.5q^2+1.1q+999, ten the reciprocal polynomial denoted A* is 
     A*=3.5+1.1q+999^2=q^2A(q^{-1})
     More formal:
     
     A^*(q)=1+a_1*q+\cdots+a_{nq}^{na}=q^{na}A(q^{-1})

     

See Åström & Wittenmark: "Computer Controlled Systems", Prentice-Hall, 1997, for a full description.

eval

public double eval(double r)

Calculates the value of the polynomial for q = r, i.e. the scalar A(r).

makeMonic

public Polynomial makeMonic()

Makes the polynomial monic by dividing all coefficients with the highest power coefficient. In a monic polynomial the highest power coefficient is always 1.

     For example let A=2q^2+q+10, then the following operation

     A.MakeMonic();

     sets A=q^2+0.5q+5.

cut

public void cut(double reps)

Normalizes the polynomial, i.e. all coefficients that are smaller than the largest coefficient times epsilon are removed. If necessary the degree of the polynomial is changed.

assign

public Polynomial assign(Polynomial p)

Gives the polynial the value of p.

shiftForward

public void shiftForward(double newLastcoeff)

This method shifts the coefficients of the polynomial to the left and sets the last coefficient to r. Let A=2.3q^2+3.5q+1.1, then the following operation

     A.ShiftForward(999);
     

changes A to A=3.5q^2+1.1q+999

shiftBackward

public void shiftBackward(double newFirstcoeff)

This method shifts the coefficients of the polynomial to the right and sets the highest power coefficient to r.\\

mul

public Polynomial mul(double r)

Returns a Polynomial whose value is (this * r).

mul

public static Polynomial mul(Polynomial p,
                             double r)

Returns a Polynomial whose value is (p * r).

mul

public static Polynomial mul(double r,
                             Polynomial p)

Returns a Polynomial whose value is (r * p).

mul

public Polynomial mul(Polynomial p2)

Returns a Polynomial whose value is (this * p2).

mul

public static Polynomial mul(Polynomial p1,
                             Polynomial p2)

Returns a Polynomial whose value is (p1 * p2).

div

public Polynomial div(double r)

Returns a Polynomial whose value is the quotient of (this / r).

div

public static Polynomial div(Polynomial p,
                             double r)

Returns a Polynomial whose value is (p / r).

div

public Polynomial div(Polynomial p2)

Returns a Polynomial whose value is the quotient of (this / p2).

div

public static Polynomial div(Polynomial p1,
                             Polynomial p2)

Returns a Polynomial whose value is the quotient of (p1 / p2).

mod

public Polynomial mod(Polynomial p2)

Returns a Polynomial whose value is the remainder of (this / p2).

mod

public static Polynomial mod(Polynomial p1,
                             Polynomial p2)

Returns a Polynomial whose value is the remainder of (p1 / p2).

add

public Polynomial add(Polynomial p2)

Returns a Polynomial whose value is the quotient of (this + p2).

add

public static Polynomial add(Polynomial p1,
                             Polynomial p2)

Returns a Polynomial whose value is the quotient of (p1 + p2).

sub

public Polynomial sub(Polynomial p2)

Returns a Polynomial whose value is the quotient of (this - p2).

sub

public static Polynomial sub(Polynomial p1,
                             Polynomial p2)

Returns a Polynomial whose value is the quotient of (p1 + p2).

isStable

public boolean isStable()

This method uses Jury's Criterion to determine if the polynomial is stable. The algorithm is described in Åström & Wittenmark: "Computer Controlled Systems", Prentice-Hall. If the polynomial is stable true is returned otherwise false.

equals

public boolean equals(Polynomial p2)

Compares this Polynomial with the specified Polynomial for equality. Returns true if and only if the Polynomials are numerically equal to each other.

GCD

public static void GCD(Polynomial a,
                       Polynomial b,
                       double reps,
                       Polynomial g,
                       Polynomial x,
                       Polynomial y,
                       Polynomial u,
                       Polynomial v)

This function calculates the greatest common divisor of polynomials A and B using Euclid's algorithm.

      a*x + b*y = g

      a*u + b*v = 0
      

g is the greatest common divisor of a and b.

diophantineMDS

public static int diophantineMDS(Polynomial a,
                                 Polynomial b,
                                 Polynomial c,
                                 Polynomial xx,
                                 Polynomial yy)

DiophantineMDS finds the minimal degree solution to the equation AR+BS=Ac. The equation is solved in two steps. First the two equations

     AX+BY=G 
     AU+BV=0 
     

are solved using the GCD function. Then in the second step the general solution is found through

     R = R0 + Q*U
     S = S0 + Q*V

     with
     
     R0 = X*Ac div G 
     S0 = Y*Ac div G
     

where G is the greatest common divisor of A and B. The minimal degree solution is now given by simply choosing Q=S0 div V. The greatest common G must divide Ac otherwise the equation has no solution.

Example

     Polynomial A, B, C, XX, YY;

     A = new Polynomial();
     B = new Polynomial();
     C = new Polynomial();
     XX = new Polynomial();
     YY = new Polynomial();

     real A_data[]= {1, 2.3, 3.5};
     A.Set(2, A_data, "A");
     
     real B_data[] = {1, 1.45};
     B.Set(1, B_data, "B");
     
     real C_data[] = {1, 3.4, 0.8, 2};
     C.Set(3, C_data, "C");
     
     DiophantineMDS(A, B, C, XX, YY);
     
     XX.Display("XX = ");
     YY.Display("YY = ");

     

When running the program above the following output is produced:

     
     XX = 1.000 q + 3.629
     YY = -2.529 q -7.379
     
     

integralDesign

public static void integralDesign(Polynomial A,
                                  Polynomial B,
                                  Polynomial Ac,
                                  Polynomial R,
                                  Polynomial S,
                                  double x0)

robustDesign

public static void robustDesign(Polynomial A,
                                Polynomial B,
                                Polynomial Ac,
                                Polynomial R,
                                Polynomial S,
                                double x0)

robustIntegralDesign

public static void robustIntegralDesign(Polynomial A,
                                        Polynomial B,
                                        Polynomial Ac,
                                        Polynomial R,
                                        Polynomial S,
                                        double x1,
                                        double x0)

LQGDesign

public static void LQGDesign(Polynomial A,
                             Polynomial B,
                             Polynomial C,
                             Polynomial R,
                             Polynomial S,
                             double rho)

Calculates a LQG-controller for the system $A$, $B$, and $C$ with the loss function coefficient $\rho$. The computational procedure is described in detail in section 12.5 in \cite{ast+wit90}. This implementation only handles the case where $A(0)\neq0$. Below is an example of how a LQG-problem can be solved using this function. It is taken from \cite{ast+wit90}, Example 12.7, with $A(z)=z+0.5$, $B(z)=0.5$, and $C=z+0.8$.

Example  
      
   Polynomial A, B, C, R, S;
   A = new Polynomial();
   B = new Polynomial();
   C = new Polynomial();
   R = new Polynomial();
   S = new Polynomial();

   System.out.println("Example 12.7 in CCS");
   System.out.println("LQG Design with rho = 0.5.");
 
   A.SetDegree(1);
   A[0] = 1; 
   A[1] = 0.5;  // a = 0.5
   B.SetDegree(0);
   B[0] = 0.5;  // b = 0.5
   C.SetDegree(1);   
   C[0] = 1; 
   C[1] =0.8;   // c = 0.8
  
   LQGDesign(A, B, C, R, S, 0.5);  
   
   R.Display("R = ");
   S.Display("S = "); 

   The following output is generated:

   R = 1.000 q + 0.614
   S = -0.112
   

MDPPNZCDesign

public static int MDPPNZCDesign(Polynomial A,
                                Polynomial B,
                                Polynomial Am,
                                Polynomial Ao,
                                Polynomial R,
                                Polynomial S,
                                Polynomial T)

MDPPNZCDesign stands for Minimal Degree Pole Placement with No Zero Cancelation. The function chooses B+ = 1 and B- = B. Furthermore Bm is chosen so that the stationary gain will be unity.

     Bm(q) = A_m(1)/B(q) B(1)
     

The closed-loop characteristic equation to be solved now becomes

     A R + B S = A_o A_m
     

The T polynomial is given by

     T(q) = A_m(1) / A_o(q) B(1)
     

The usage of MDPPNZCDesign-method is used is demonstrated below.

 Example
     
     Polynomial A, B, Am, Ao, R, S, T;
     A = new Polynomial();
     B = new Polynomial();
     Am = new Polynomial();
     Ao = new Polynomial();
     R = new Polynomial();
     S = new Polynomial();
     T = new Polynomial();

     B.SetDegree(1);
     B[0] = 1;
     B[1] = 0.7;
     A.SetDegree(2);
     A[0] = 1;
     A[1] = -1.8;
     A[2] = 0.81;
     Am.SetDegree(2);
     Am[0] = 1;
     Am[1] = -1.5;
     Am[2] = 0.7;
     Ao.SetDegree(1);
     Ao[0] = 1;
     Ao[1] = 0;
     
     MDPPNZCDesign(A, B, Am, Ao, R, S, T);
     
     R.Display("R = ");
     S.Display("S = ");
     T.Display("T = ");
     }

     From this program the following output is generated:
     
     R = 1.000 q + 0.088
     S = 0.213 q + -0.101
     T = 0.118 q + 0.000
     

sfactorize

public static void sfactorize(Polynomial b,
                              Polynomial a)

The reciprocal polynomial of A, denoted A* is obtained by reversing the order of the coefficients of A.

 Example
     Let  A(z) = z^(na)+a_1*z^(na-1) + ... + a_(na).
     
     The it reciprocal polynomial of A is the given as
     
     A*(z) = 1 + a_1*z + ...  + a_(na) * z^(na) = z^(na) * A(z^(-1)).
     

Sfactorize takes a polynomial B = e e* and returns a stable polynomial A so that A A* = e e*

dyadicReduction

public static void dyadicReduction(double[] a,
                                   double[] b,
                                   double[] alpha,
                                   double[] beta,
                                   int i0,
                                   int i1,
                                   int i2)

Given vectors

     a = [1 a_2 ... a_n]^T
     b = [1 b_2 ... b_n]^T
     

and scalars alpha and beta, find vectors

     a' = [1 a'_2 ... a'_n]^T
     b' = [1 b'_2 ... b'_n]^T
     

such that

     alpha a a^T + \beta b b^T = alpha' a' a'^T + beta' b' b'^T
     

The vectors a' and b' can be found using dyadic decomposition. DyadicReduction is very useful when doing square root recursive least square estimations. For a closer look at the algorithm and its applications see Åström & Wittenmark: "Adaptive Control", Addison-Wesley, 1995.

LDFilter

public static void LDFilter(double[] l,
                            double[] d,
                            double[] phi,
                            double[] theta,
                            double lambda,
                            int degree)

LFilter s an implementation of an recursive estimator and would together with any of the design functions above form a complete adaptive controller. The algorithms behind LDFilter is taken from Åström & Wittenmark: "Adaptive Control", Addison-Wesley, 1995. Let the system to be estimated be on the form

  y(t) = phi^T(t) theta 
  

where theta is a parameter vector and phi is a vector of signals. The following recursive least-square estimation algorithm is used

  theta'(t) = theta'(t-1) + K(t)(y(t)-phi^T(t)theta'(t-1)
  K(t) = P(t) phi(t)
  P(t) =(I-K(t) phi^T(t))P(t-1)
  

The covariance matrix P has a decomposition P=LDL^T, where L is a lower triangular matrix and D is a diagonal matrix. Initially set L=I, which gives that P(0)=D.

Parameters:

l - A matrix of doubles with the size (degree+1)^2.

d - An array of doubles with the size (degree+1).

Phi - The Phi vector contains old process values and old control signals. The polynomial is arranged on the following format:

  Phi = [y(t), -y(t-1), \ldots , -y(t-n), u(t-d), \dots , u(t-d-m)]
  

where n=deg(A), m=deg(B) and d=n-m and is of size (degree+1).

theta - This is a vector with the size set to (deg(A) + deg(B) + 2) and with the coefficients to be estimated on the following format:

  theta = [a_0, a_1, \ldots, a_n,b_0, ..., b_m]
  

and is of size (degree+1).

l - This parameter is the forgetting factor $\lambda$.

degree - This parameter is the dimension of the problem degree = deg(A) + deg(B) + 1.

set

public void set(int newDeg,
                double[] newCoeff)

Sets the degree of the Polynomial to newDeg and uses the values in newCoeff as coefficients.

getDegree

public int getDegree()

Returns the degree of the polynomial.

get

public double get(int elem)

Returns the value of element i. A polynomial

     A=a_0 z^n + a_1 z^(n-1) + ... + a_n

     is organized as  
     A=[a_0, a_1, ... ,a_n]
     

This means that for example A.get(0) will return the value of the highest power coefficient.

set

public void set(int elem,
                double value)

Sets the value of element i. A polynomial

     A=a_0 z^n + a_1 z^(n-1) + ... + a_n

     is organized as  
     A=[a_0, a_1, ... ,a_n]
     

This means that for example A.set(0, 1) will set the value of the highest power coefficient to 1.