diff options
Diffstat (limited to '2048')
36 files changed, 848 insertions, 0 deletions
diff --git a/2048/DEPENDENCIES/.png b/2048/DEPENDENCIES/.png Binary files differnew file mode 100644 index 000000000..153470f38 --- /dev/null +++ b/2048/DEPENDENCIES/.png diff --git a/2048/DEPENDENCIES/cl.sci b/2048/DEPENDENCIES/cl.sci new file mode 100644 index 000000000..561874cb8 --- /dev/null +++ b/2048/DEPENDENCIES/cl.sci @@ -0,0 +1,28 @@ +// Calculation of closed loop transfer functions
+// 11.6
+
+// function [Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar] = ...
+// cl(A,dA,B,dB,C,dC,k,S,dS,R,dR,int)
+// int>=1 means integrated noise and control law:
+// delta u = - (S/R)y
+// Evaluates the closed loop transfer function and
+// variances of input and output
+
+function [Nu,dNu,Du,dDu,Ny,dNy,Dy,dDy,yvar,uvar] = ...
+ cl(A,dA,B,dB,C,dC,k,S,dS,R,dR,int1)
+[zk,dzk] = zpowk(k);
+
+[BS,dBS] = polmul(B,dB,S,dS);
+[zBS,dzBS] = polmul(zk,dzk,BS,dBS);
+[RA,dRA] = polmul(R,dR,A,dA);
+if int1>=1, [RA,dRA] = polmul(RA,dRA,[1 -1],1); end
+
+[D,dD] = poladd(RA,dRA,zBS,dzBS);
+
+[Ny,dNy] = polmul(C,dC,R,dR);
+[Nu,dNu] = polmul(C,dC,S,dS);
+
+[Nu,dNu,Du,dDu,uvar] = tfvar(Nu,dNu,D,dD);
+[Ny,dNy,Dy,dDy,yvar] = tfvar(Ny,dNy,D,dD);
+
+endfunction;
diff --git a/2048/DEPENDENCIES/desired.sci b/2048/DEPENDENCIES/desired.sci new file mode 100644 index 000000000..d84e6c918 --- /dev/null +++ b/2048/DEPENDENCIES/desired.sci @@ -0,0 +1,11 @@ +// Calculation of desired closed loop characteristic polynomial, as discussed in Sec. 7.7.
+// 9.4
+
+// function [phi,dphi] = desired(Ts,rise,epsilon)
+// Based on transient requirements,
+// calculates closed loop characteristic polynomial
+
+function [phi,dphi] = desired(Ts,rise,epsilon)
+Nr = rise/Ts; omega = %pi/2/Nr; rho = epsilon^(omega/%pi);
+phi = [1 -2*rho*cos(omega) rho^2]; dphi = length(phi)-1;
+endfunction;
diff --git a/2048/DEPENDENCIES/filtval.sci b/2048/DEPENDENCIES/filtval.sci new file mode 100644 index 000000000..1d4dcb6b5 --- /dev/null +++ b/2048/DEPENDENCIES/filtval.sci @@ -0,0 +1,11 @@ +// Value of polynomial p(x), evaluated at x
+// 11.16
+
+// finds the value of a polynomial in powers of z^{-1}
+// function Y = filtval(P,z)
+
+function Y = filtval(P,z)
+N = length(P)-1;
+Q = polyno(P,'x');
+Y = horner(Q,z)/z^N;
+endfunction;
diff --git a/2048/DEPENDENCIES/flip.sci b/2048/DEPENDENCIES/flip.sci new file mode 100644 index 000000000..c19564c9a --- /dev/null +++ b/2048/DEPENDENCIES/flip.sci @@ -0,0 +1,4 @@ +// 10.5
+function b = flip(a)
+b = a(length(a):-1:1);
+endfunction;
diff --git a/2048/DEPENDENCIES/gmv.sci b/2048/DEPENDENCIES/gmv.sci new file mode 100644 index 000000000..1e3ee068c --- /dev/null +++ b/2048/DEPENDENCIES/gmv.sci @@ -0,0 +1,17 @@ +// General minimum variance controller design, as given by Eq. 11.66 on page 421 and Eq. 11.70 on page 422.
+// 11.12
+
+// function [Sc,dSc,Rc,dRc] = gmv(A,dA,B,dB,C,dC,k,rho,int)
+// implements the generalized minimum variance controller
+// if int>=1, integrated noise is assumed; otherwise,
+// it is not integrated noise
+
+function [Sc,dSc,R,dR] = gmv(A,dA,B,dB,C,dC,k,rho,int1)
+zk = zeros(1,k+1); zk(k+1) = 1;
+if int1>=1, [A,dA] = polmul([1 -1],1,A,dA); end
+[Fk,dFk,Ek,dEk] = xdync(zk,k,A,dA,C,dC);
+[Gk,dGk] = polmul(Ek,dEk,B,dB);
+alpha0 = Gk(1)/C(1);
+Sc = alpha0 * Fk; dSc = dFk;
+[R,dR] = poladd(alpha0*Gk,dGk,rho*C,dC);
+endfunction;
diff --git a/2048/DEPENDENCIES/gmv_mac1.dat b/2048/DEPENDENCIES/gmv_mac1.dat Binary files differnew file mode 100644 index 000000000..42d664411 --- /dev/null +++ b/2048/DEPENDENCIES/gmv_mac1.dat diff --git a/2048/DEPENDENCIES/gmvc_pid.sci b/2048/DEPENDENCIES/gmvc_pid.sci new file mode 100644 index 000000000..c00f2b000 --- /dev/null +++ b/2048/DEPENDENCIES/gmvc_pid.sci @@ -0,0 +1,42 @@ +// PID tuning through GMVC law
+// 11.17
+
+// function [Kc,tau_i,tau_d,L] = gmvc_pid(A,B,k,T,Ts)
+// Determines p,i,d tuning parameters using GMVC
+// Plant model: Integrated white noise
+// A, B in discrete time form
+
+function [Kc,tau_i,tau_d,L] = gmvc_pid(A,B,k,T,Ts)
+
+dA = length(A)-1; dB = length(B)-1;
+dT = length(T)-1;
+if dA > 2,
+ disp('degree of A cannot be more than 2')
+ exit
+elseif dB > 1,
+ disp('degree of B cannot be more than 1')
+ exit
+elseif dT > 2,
+ disp('degree of T cannot be more than 2')
+ exit
+end
+delta = [1 -1]; ddelta = 1;
+
+[Adelta,dAdelta] = polmul(A,dA,delta,ddelta);
+
+[Q,dQ,P,dP] = ...
+xdync(Adelta,dAdelta,B,dB,T,dT);
+PAdelta = P(1)*Adelta;
+
+[zk,dzk] = zpowk(k);
+[E,degE,F,degF] = ...
+xdync(PAdelta,dAdelta,zk,dzk,P,dP);
+nu = P(1)*E(1)*B(1);
+Kc = -1/nu*(F(2)+2*F(3));
+tau_i = -(F(2)+2*F(3))/(F(1)+F(2)+F(3))*Ts;
+tau_d = -F(3)/(F(2)+2*F(3))*Ts;
+L(1) = 1+Ts/tau_i+tau_d/Ts;
+L(2) = -(1+2*tau_d/Ts);
+L(3) = tau_d/Ts;
+L = Kc * L';
+endfunction;
diff --git a/2048/DEPENDENCIES/gpc_N.sci b/2048/DEPENDENCIES/gpc_N.sci new file mode 100644 index 000000000..2abbe0f74 --- /dev/null +++ b/2048/DEPENDENCIES/gpc_N.sci @@ -0,0 +1,31 @@ +// Calculates the GPC law given by Eq. 12.36 on page 446.
+// 12.5
+
+function [K,KH1,KH2,Tc,dTc,Sc,dSc,R1,dR1] = ...
+gpc_N(A,dA,B,dB,k,N1,N2,Nu,rho)
+D=[1 -1]; dD=1; AD=convol(A,D); dAD=dA+1;
+zj = 1; dzj = 0;
+for i = 1:N1+k-1
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+end
+G = zeros(N2-N1+1,Nu+1);
+H1 = zeros(N2-N1+1,k-1+dB); H2 = zeros(N2-N1+1,dA+1);
+for j = k+N1:k+N2
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+ [Fj,dFj,Ej,dEj] = xdync(zj,dzj,AD,dAD,1,0);
+ [Gj,dGj] = polmul(B,dB,Ej,dEj);
+ if (j-k >= Nu)
+ G(j-(k+N1-1),1:Nu+1) = flip(Gj(j-k-Nu+1:j-k+1));
+else
+ G(j-(k+N1-1),1:j-k+1) = flip(Gj(1:j-k+1));
+end
+ H1(j-(k+N1-1),1:k-1+dB) = Gj(j-k+2:j+dB);
+ H2(j-(k+N1-1),1:dA+1) = Fj;
+end
+K = inv(G'*G+rho*eye(Nu+1,Nu+1))*G';
+// Note: inverse need not be calculated
+KH1 = K * H1; KH2 = K * H2;
+R1 = [1 KH1(1,:)]; dR1 = length(R1)-1;
+Sc = KH2(1,:); dSc = length(Sc)-1;
+Tc = K(1,:); dTc = length(Tc)-1;
+endfunction;
diff --git a/2048/DEPENDENCIES/gpc_Nc.sci b/2048/DEPENDENCIES/gpc_Nc.sci new file mode 100644 index 000000000..9401ffde1 --- /dev/null +++ b/2048/DEPENDENCIES/gpc_Nc.sci @@ -0,0 +1,37 @@ +// Calculates the GPC law for different prediction and control horizons
+// 12.9
+
+function [K,KH1,KH2,Tc,dTc,Sc,dSc,R1,dR1] = ...
+gpc_Nc(A,dA,B,dB,C,dC,k,N1,N2,Nu,rho)
+D=[1 -1]; dD=1; AD=convol(A,D); dAD=dA+1;
+zj = 1; dzj = 0;
+for i = 1:N1+k-1
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+end
+M = 2*k+N2-2+dB; P = max(k+N2+dA-1,dC-1)
+G = zeros(N2-N1+1,Nu+1); H1 = zeros(N2-N1+1,M);
+H2 = zeros(N2-N1+1,P+1);
+for j = k+N1:k+N2
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+ [Fj,dFj,Ej,dEj] = xdync(zj,dzj,AD,dAD,C,dC);
+ [Nj,dNj,Mj,dMj] = xdync(zj,dzj,C,dC,1,0);
+ [Gj,dGj] = polmul(Mj,dMj,Ej,dEj);
+ [Gj,dGj] = polmul(Gj,dGj,B,dB);
+ [Pj,dPj] = polmul(Mj,dMj,Fj,dFj);
+ [Pj,dPj] = poladd(Nj,dNj,Pj,dPj);
+ if (j-k >= Nu)
+ G(j-(k+N1-1),1:Nu+1) = flip(Gj(j-k-Nu+1:j-k+1));
+else
+ G(j-(k+N1-1),1:j-k+1) = flip(Gj(1:j-k+1));
+end
+ H1(j-(k+N1-1),1:j+k-2+dB) = Gj(j-k+2:2*j+dB-1);
+ dPj = max(j-1+dA,dC-1);
+ H2(j-(k+N1-1),1:dPj+1) = Pj;
+end
+K = inv(G'*G+rho*eye(Nu+1,Nu+1))*G';
+// Note: inverse need not be calculated
+KH1 = K * H1; KH2 = K * H2;
+R1 = [1 KH1(1,:)]; dR1 = length(R1)-1;
+Sc = KH2(1,:); dSc = length(Sc)-1;
+Tc = K(1,:); dTc = length(Tc)-1;
+endfunction;
diff --git a/2048/DEPENDENCIES/gpc_bas.sci b/2048/DEPENDENCIES/gpc_bas.sci new file mode 100644 index 000000000..dd742e0c5 --- /dev/null +++ b/2048/DEPENDENCIES/gpc_bas.sci @@ -0,0 +1,23 @@ +// Calculates the GPC law given by Eq. 12.19 on page 441.
+// 12.2
+
+function [K,KH1,KH2,Tc,dTc,Sc,dSc,R1,dR1] = ...
+gpc_bas(A,dA,B,dB,N,k,rho)
+D=[1 -1]; dD=1; AD=convol(A,D); dAD=dA+1; Nu=N+1;
+zj = 1; dzj = 0; G = zeros(Nu,Nu);
+H1 = zeros(Nu,k-1+dB); H2 = zeros(Nu,dA+1);
+for j = 1:Nu,
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+ [Fj,dFj,Ej,dEj] = xdync(zj,dzj,AD,dAD,1,0);
+ [Gj,dGj] = polmul(B,dB,Ej,dEj);
+ G(j,1:dGj) = flip(Gj(1:dGj));
+ H1(j,1:k-1+dB) = Gj(dGj+1:dGj+k-1+dB);
+ H2(j,1:dA+1) = Fj;
+end
+K = inv(G'*G+rho*eye(Nu,Nu))*G';
+// Note: inverse need not be calculated
+KH1 = K * H1; KH2 = K * H2;
+R1 = [1 KH1(1,:)]; dR1 = length(R1)-1;
+Sc = KH2(1,:); dSc = length(Sc)-1;
+Tc = K(1,:); dTc = length(Tc)-1;
+endfunction;
diff --git a/2048/DEPENDENCIES/gpc_col.sci b/2048/DEPENDENCIES/gpc_col.sci new file mode 100644 index 000000000..569718f08 --- /dev/null +++ b/2048/DEPENDENCIES/gpc_col.sci @@ -0,0 +1,30 @@ +// Calculates the GPC law given by Eq. 12.36 on page 446.
+// 12.6
+
+function [K,KH1,KH2,Tc,dTc,Sc,dSc,R1,dR1] = ...
+gpc_col(A,dA,B,dB,C,dC,N,k,rho)
+D=[1 -1]; dD = 0; AD=convol(A,D); dAD=dA+1; zj=1; dzj=0;
+Nu = N+1; G=zeros(Nu,Nu); H1=zeros(Nu,2*k+N-2+dB);
+H2 = zeros(Nu,k+N+dA);
+for j = 1:Nu,
+ zj = convol(zj,[0,1]); dzj = dzj + 1;
+ [Fj,dFj,Ej,dEj] = ...
+ xdync(zj,dzj,AD,dAD,C,dC);
+ [Nj,dNj,Mj,dMj] = ...
+ xdync(zj,dzj,C,dC,1,0);
+ [Gj,dGj] = polmul(Mj,dMj,Ej,dEj);
+ [Gj,dGj] = polmul(Gj,dGj,B,dB);
+ [Pj,dPj] = polmul(Mj,dMj,Fj,dFj);
+ [Pj,dPj] = poladd(Nj,dNj,Pj,dPj);
+ j,Fj,Ej,Mj,Nj,Gj,Pj
+ G(j,1:j) = flip(Gj(1:j));
+ H1(j,1:dGj-j+1) = Gj(j+1:dGj+1);
+ H2(j,1:dPj+1) = Pj;
+end
+K = inv(G'*G+rho*eye(Nu,Nu))*G'
+// Note: inverse need not be calculated
+KH1 = K * H1; KH2 = K * H2;
+R1 = [1 KH1(1,:)]; dR1 = length(R1)-1;
+Sc = KH2(1,:); dSc = length(Sc)-1;
+Tc = K(1,:); dTc = length(Tc)-1;
+endfunction;
diff --git a/2048/DEPENDENCIES/gpc_pid.sci b/2048/DEPENDENCIES/gpc_pid.sci new file mode 100644 index 000000000..2359a5dfe --- /dev/null +++ b/2048/DEPENDENCIES/gpc_pid.sci @@ -0,0 +1,52 @@ +// Predictive PID, tuned with GPC, as explained in Sec. 12.2.3.
+// 12.11
+
+function [Kp,Ki,Kd] = ...
+gpc_pid(A,dA,B,dB,C,dC,N1,N2,Nu,lambda,gamm,gamma_y)
+Adelta=convol(A,[1 -1]); G=[];
+for i=N1:N2
+ zi=zpowk(i);
+ [E,dE,F,dF]=xdync(Adelta,dA+1,zi,i,C,dC);
+ [Gtilda,dGtilda,Gbar,dGbar] = ...
+ xdync(C,dC,zi,i,E*B,dE+dB);
+ for j = 1:i, Gtilda1(j)=Gtilda(i+1-j); end
+ Gtilda2 = Gtilda1.'; // Added because Scilab forms a column vecor
+ // while Matlab forms a row vector, by default
+ if i<=Nu-1
+ G=[G;[Gtilda2,zeros(1,Nu-i)]];
+ else
+ G=[G;Gtilda2(1:Nu)];
+ end
+end
+es=sum(C)/sum(A); gs=sum(B)/sum(A); F_s=es*A; G_s=[];
+for i=1:Nu
+ if ((Nu - i) == 0)
+ row=gs*ones(1,i);
+ else
+ row=gs*ones(1,i); row=[row,zeros(Nu-i,Nu-i)];
+ end;
+ G_s=[G_s;row];
+end
+lambda_mat=lambda*(diag(ones(1,Nu)));
+gamma_mat=gamm*(diag(ones(1,Nu)));
+gamma_y_mat=gamma_y*(diag(ones(1,N2-N1+1)));
+mat1=inv(G'*gamma_y_mat*G+lambda_mat+G_s'*gamma_mat*G_s);
+mat2=mat1*(G'*gamma_y_mat);
+mat2_s=mat1*(G_s'*gamma_mat);
+h_s=sum(mat2_s(1,:)); h=mat2(1,:);
+T=C; R=C*(sum(h(:))+h_s); S=0;
+for i=N1:N2
+ zi=zpowk(i);
+ [E,dE,F,dF]=xdync(Adelta,dA+1,zi,i,C,dC);
+ [Gtilda,dGtilda,Gbar,dGbar]=...
+ xdync(C,dC,zi,i,E*B,dE+dB);
+ S=S+F*h(i);
+end
+S=S+F_s*h_s;
+if length(A)==3
+ Kp=S(1)-R-S(3); Ki=R; Kd=S(3);
+else
+ Kp=S(1)-R; Ki=R; Kd=0;
+end
+
+endfunction;
diff --git a/2048/DEPENDENCIES/imc_stable.sci b/2048/DEPENDENCIES/imc_stable.sci new file mode 100644 index 000000000..d53da9c82 --- /dev/null +++ b/2048/DEPENDENCIES/imc_stable.sci @@ -0,0 +1,31 @@ +// Design of conventional controller which is an equivalent of internal model controller
+// 10.9
+
+// Designs Discrete Internal Model Controller
+// for transfer function z^{-k}B(z^{-1})/A(z^{-1})
+// Numerator and Denominator of IMC HQ are outputs
+// Controller is also given in R,S form
+
+
+function [k,HiN,HiD,R,S,mu] = imc_stable(B,A,k,alpha)
+
+[Kp,d,Bg,Bnmp,Bm] = imcsplit(B,mtlb_logical(1));
+Bg = Kp * Bg;
+
+Bnmpr = flip(Bnmp);
+Bms = sum(Bm);
+HiN = A;
+HiD = Bms * convol(Bg,Bnmpr);
+k = k+d;
+
+[zk,dzk] = zpowk(k);
+Bf = (1-alpha);
+Af = [1 -alpha];
+S = convol(Bf,A);
+R1 = convol(Af,convol(Bnmpr,Bms));
+R2 = convol(zk,convol(Bf,convol(Bnmp,Bm)));
+
+[R,dR] = poladd(R1,length(R1)-1,-R2,length(R2)-1);
+R = convol(Bg,R);
+endfunction;
+
diff --git a/2048/DEPENDENCIES/imc_stable1.sci b/2048/DEPENDENCIES/imc_stable1.sci new file mode 100644 index 000000000..6b8dec788 --- /dev/null +++ b/2048/DEPENDENCIES/imc_stable1.sci @@ -0,0 +1,19 @@ +// Design of internal model controller
+// 10.4
+// Designs Discrete Internal Model Controller
+// for transfer function z^{-k}B(z^{-1})/A(z^{-1})
+// Numerator and Denominator of IMC HQ are outputs
+// Controller is also given in R,S form
+
+function [k,HiN,HiD] = imc_stable1(B,A,k,alpha)
+
+[Kp,d,Bg,Bnmp,Bm] = imcsplit(B,mtlb_logical(1));
+Bg = Kp * Bg;
+Bnmpr = flip(Bnmp);
+Bms = sum(Bm);
+HiN = A;
+HiD = Bms * convol(Bg,Bnmpr);
+k = k+d;
+endfunction;
+
+
diff --git a/2048/DEPENDENCIES/imcsplit.sci b/2048/DEPENDENCIES/imcsplit.sci new file mode 100644 index 000000000..59e6ef95f --- /dev/null +++ b/2048/DEPENDENCIES/imcsplit.sci @@ -0,0 +1,37 @@ +// Splitting a polynomial B(z)
+// 10.3
+// Splits a polynomial B into good, nonminimum with
+// positive real & with negative real parts.
+// All are returned in polynomial form.
+// Gain is returned in Kp and delay in k.
+
+function [Kp,k,Bg,Bnmp,Bm] = imcsplit(B,polynomial)
+k = 0;
+Kp = 1;
+if(polynomial)
+ rts = roots(B);
+ Kp = sum(B)/sum(coeff(poly(rts,'z')));
+else
+ rts = B;
+end
+Bg = 1; Bnmp = 1; Bm = 1;
+for i = 1:length(rts),
+ rt = rts(i);
+ if rt == 0,
+ k = k+1;
+ elseif (abs(rt)<1 & real(rt)>=0)
+ Bg = convol(Bg,[1 -rt]);
+ elseif (abs(rt)>=1 & real(rt)>=0)
+ Bnmp = convol(Bnmp,[1 -rt]);
+ else
+ Bm = convol(Bm,[1 -rt]);
+ end
+end
+
+
+
+
+
+
+
+
diff --git a/2048/DEPENDENCIES/kal_ex.sci b/2048/DEPENDENCIES/kal_ex.sci new file mode 100644 index 000000000..739e96537 --- /dev/null +++ b/2048/DEPENDENCIES/kal_ex.sci @@ -0,0 +1,12 @@ +// Kalman filter example of estimating a constant
+// 14.4
+
+function [xhat,P,y] = kal_ex(x,xline,M)
+y = x + rand();
+Q = 0; R = 1;
+xhat_ = xline;
+P_ = M + Q;
+K = P_/(P_+R);
+P = (1-K)*P_;
+xhat = xhat_ + K*(y-xhat_);
+endfunction;
diff --git a/2048/DEPENDENCIES/lqg1.sci b/2048/DEPENDENCIES/lqg1.sci new file mode 100644 index 000000000..2ce4a74ec --- /dev/null +++ b/2048/DEPENDENCIES/lqg1.sci @@ -0,0 +1,48 @@ +// LQG control design by polynomial method, to solve Eq. 13.51 on page 472.
+// 13.4
+
+// LQG controller design by method of Ahlen and Sternad
+// function [R,degR,S,degS] = ...
+// lqg(A,degA,B,degB,C,degC,k,rho,V,degV,W,degW,F,degF)
+
+function [R,degR,S,degS] = ...
+lqg1(A,degA,B,degB,C,degC,k,rho,V,degV,W,degW,F,degF)
+
+[r,b,degb] = ...
+specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF);
+
+WFA = flip(convol(A,convol(F,W)));
+dWFA = degW + degF + degA;
+
+[rhs1,drhs1] = polmul(W,degW,WFA,dWFA);
+[rhs1,drhs1] = polmul(rhs1,drhs1,C,degC);
+rhs1 = rho * rhs1;
+rhs2 = convol(C,convol(V,flip(convol(B,V))));
+drhs2 = degC + 2*degV + degB;
+for i = 1:degb-degB-degV,
+ rhs2 = convol(rhs2,[0,1]);
+end
+drhs2 = drhs2 + degb-degB-degV;
+C1 = zeros(1,2);
+
+[C1,degC1] = putin(C1,0,rhs1,drhs1,1,1);
+[C1,degC1] = putin(C1,degC1,rhs2,drhs2,1,2);
+rbf = r * flip(b);
+D1 = zeros(2,2);
+[D1,degD1] = putin(D1,0,rbf,degb,1,1);
+for i = 1:k,
+ rbf = convol(rbf,[0 1]);
+end
+[D1,degD1] = putin(D1,degD1,rbf,degb+k,2,2);
+N = zeros(1,2);
+[N,degN] = putin(N,0,-B,degB,1,1);
+[AF,dAF] = polmul(A,degA,F,degF);
+[N,degN] = putin(N,degN,AF,dAF,1,2);
+
+[Y,degY,X,degX] = xdync(N,degN,D1,degD1,C1,degC1);
+
+[R,degR] = ext(X,degX,1,1);
+[S,degS] = ext(X,degX,1,2);
+X = flip(Y);
+
+endfunction;
diff --git a/2048/DEPENDENCIES/lqg_visc.dat b/2048/DEPENDENCIES/lqg_visc.dat Binary files differnew file mode 100644 index 000000000..888d22c6b --- /dev/null +++ b/2048/DEPENDENCIES/lqg_visc.dat diff --git a/2048/DEPENDENCIES/mv.sci b/2048/DEPENDENCIES/mv.sci new file mode 100644 index 000000000..532fc6b4d --- /dev/null +++ b/2048/DEPENDENCIES/mv.sci @@ -0,0 +1,16 @@ +// Minimum variance control law design, given by Eq. 11.40 on page 413.
+// 11.5
+
+// function [S,dS,R,dR] = mv(A,dA,B,dB,C,dC,k,int)
+// implements the minimum variance controller
+// if int>=1, integrated noise is assumed; otherwise,
+// it is not integrated noise
+
+function [S,dS,R,dR] = mv(A,dA,B,dB,C,dC,k,int1)
+zk = zeros(1,k+1); zk(k+1) = 1;
+if int1>=1, [A,dA] = polmul([1 -1],1,A,dA); end
+[Fk,dFk,Ek,dEk] = xdync(zk,k,A,dA,C,dC);
+
+[Gk,dGk] = polmul(Ek,dEk,B,dB);
+S = Fk; dS = dFk; R = Gk; dR = dGk;
+endfunction;
diff --git a/2048/DEPENDENCIES/mv_nm.sci b/2048/DEPENDENCIES/mv_nm.sci new file mode 100644 index 000000000..23b902f45 --- /dev/null +++ b/2048/DEPENDENCIES/mv_nm.sci @@ -0,0 +1,16 @@ +// Minimum variance control for nonminimum phase systems
+// 11.9
+
+// function [Sc,dSc,Rc,dRc] = mv_mv(A,dA,B,dB,C,dC,k,int)
+// implements the minimum variance controller
+// if int>=1, integrated noise is assumed; otherwise,
+// it is not integrated noise
+
+function [Sc,dSc,Rc,dRc] = mv_nm(A,dA,B,dB,C,dC,k,int1)
+if int1>=1, [A,dA] = polmul([1 -1],1,A,dA); end
+[zk,dzk] = zpowk(k);
+[Bzk,dBzk] = polmul(B,dB,zk,dzk);
+[Bg,Bb] = polsplit3(B); Bbr = flip(Bb);
+RHS = convol(C,convol(Bg,Bbr)); dRHS = length(RHS)-1;
+[Sc,dSc,Rc,dRc] = xdync(Bzk,dBzk,A,dA,RHS,dRHS);
+endfunction;
diff --git a/2048/DEPENDENCIES/myc2d.sci b/2048/DEPENDENCIES/myc2d.sci new file mode 100644 index 000000000..1644cddfe --- /dev/null +++ b/2048/DEPENDENCIES/myc2d.sci @@ -0,0 +1,20 @@ +// Discretization of continuous transfer function. The result is numerator and denominator in powers of z^{-1} and the delay term k.
+// 9.2
+// function [B,A,k] = myc2d(G,Ts)
+// Produces numerator and denominator of discrete transfer
+// function in powers of z^{-1}
+// G is continuous transfer function; time delays are not allowed
+// Ts is the sampling time, all in consistent time units
+
+function [B,A,k] = myc2d(G,Ts)
+H = ss2tf(dscr(G,Ts));
+num1 = coeff(H('num'));
+den1 = coeff(H('den'));//-------------
+A = den1(length(den1):-1:1);
+num2 = num1(length(num1):-1:1); //flip
+nonzero = find(num1);
+first_nz = nonzero(1);
+B = num2(first_nz:length(num2)); //-------------
+k = length(den1) - length(num1);
+endfunction
+
diff --git a/2048/DEPENDENCIES/pacf.sci b/2048/DEPENDENCIES/pacf.sci new file mode 100644 index 000000000..df383297c --- /dev/null +++ b/2048/DEPENDENCIES/pacf.sci @@ -0,0 +1,24 @@ +// Determination of the PACF of AR(p) process, as explained in Sec. 6.4.5.
+// 6.10
+
+function [ajj] = pacf(v,M)
+exec('label.sci',-1);
+rvvn = xcorr(v,'coeff');
+len = length(rvvn);
+zero = (len+1)/2;
+rvvn0 = rvvn(zero);
+rvvn_one_side = rvvn(zero+1:len);
+ajj = [];
+exec('pacf_mat.sci',-1);
+for j = 1:M,
+ ajj = [ajj pacf_mat(rvvn0,rvvn_one_side,j,1)];
+end
+p = 1:length(ajj);
+N = length(p);
+lim = 2/sqrt(length(v));
+
+// Plot the figure
+plot(p,ajj,p,ajj,'o',p,lim*ones(N,1),'--',...
+ p,-lim*ones(N,1),'--');
+label('',4,'Lag','PACF',4);
+endfunction;
diff --git a/2048/DEPENDENCIES/pacf_mat.sci b/2048/DEPENDENCIES/pacf_mat.sci new file mode 100644 index 000000000..ba250b63f --- /dev/null +++ b/2048/DEPENDENCIES/pacf_mat.sci @@ -0,0 +1,44 @@ +// Construction of square matrix required to compute PACF ajj, useful for the calculations in Sec. 6.4.5.
+// 6.11
+
+function ajj = pacf_mat(rvv0,rvv_rest,p,k)
+if argn(2) == 3,
+ k = 1;
+end
+for i = 1:p
+ for j = 1:p
+ index = (k+i-1)-j;
+ if index == 0,
+ A(i,j) = rvv0;
+ elseif index < 0,
+ A(i,j) = rvv_rest(-index);
+ else
+ A(i,j) = rvv_rest(index);
+ end
+ end
+ b(i) = -rvv_rest(k+i-1);
+end
+a = A\b;
+ajj = a(p);
+endfunction;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/2048/DEPENDENCIES/pd.sci b/2048/DEPENDENCIES/pd.sci new file mode 100644 index 000000000..110401281 --- /dev/null +++ b/2048/DEPENDENCIES/pd.sci @@ -0,0 +1,14 @@ +// PD control law from polynomial coefficients, as explained in Sec. 9.8.
+// 9.22
+
+function [K,taud,N] = pd(Rc,Sc,Ts)
+
+// Both Rc and Sc have to be degree one polynomials
+
+s0 = Sc(1); s1 = Sc(2);
+r1 = Rc(2);
+K = (s0+s1)/(1+r1);
+N = (s1-s0*r1)/r1/(s0+s1);
+taudbyN = -Ts*r1/(1+r1);
+taud = taudbyN * N;
+endfunction;
diff --git a/2048/DEPENDENCIES/polsplit2.sci b/2048/DEPENDENCIES/polsplit2.sci new file mode 100644 index 000000000..524c35442 --- /dev/null +++ b/2048/DEPENDENCIES/polsplit2.sci @@ -0,0 +1,36 @@ +// Procedure to split a polynomial into good and bad factors, as discussed in Sec. 9.2.
+// 9.3
+// function [goodpoly,badpoly] = polsplit2(fac,a)
+// Splits a scalar polynomial of z^{-1} into good and bad
+// factors.
+// Input is a polynomial in increasing degree of z^{-1}
+// Optional input is a, where a <= 1.
+// Factor that has roots of z^{-1} outside a is called
+// good and the rest bad.
+// If a is not specified, it will be assumed as 1-1.0e-5
+
+function [goodpoly,badpoly] = polsplit2(fac,a)
+if argn(2) == 1, a = 1-1.0e-5; end
+if a>1 error('good polynomial is unstable'); end
+fac1 = poly(fac(length(fac):-1:1),'z','coeff');
+rts1 = roots(fac1);
+rts = rts1(length(rts1):-1:1);
+
+// extract good and bad roots
+badindex = find(abs(rts)>=a); // mtlb_find has been replaced by find
+badpoly = coeff(poly((rts(badindex)),"z","roots"));
+goodindex = find(abs(rts)<a); // mtlb_find has been replaced by find
+goodpoly = coeff(poly(rts(goodindex),"z","roots"));
+
+// scale by equating the largest terms
+[m,index] = max(abs(fac));
+goodbad = convol(goodpoly,badpoly);
+goodbad1 = goodbad(length(goodbad):-1:1);//--
+factor1 = fac(index)/goodbad1(index);//--
+goodpoly = goodpoly * factor1;
+goodpoly = goodpoly(length(goodpoly):-1:1);
+badpoly = badpoly(length(badpoly):-1:1);
+endfunction;
+
+
+
diff --git a/2048/DEPENDENCIES/polsplit3.sci b/2048/DEPENDENCIES/polsplit3.sci new file mode 100644 index 000000000..caffc7509 --- /dev/null +++ b/2048/DEPENDENCIES/polsplit3.sci @@ -0,0 +1,33 @@ +// Procedure to split a polynomial into good and bad factors, as discussed in Sec. 9.5. The factors that have roots outside unit circle or with negative real parts are defined as bad.
+// 9.12
+
+// function [goodpoly,badpoly] = polsplit3(fac,a)
+// Splits a scalar polynomial of z^{-1} into good and bad
+// factors. Input is a polynomial in increasing degree of
+// z^{-1}. Optional input is a, where a <= 1.
+// Factors that have roots outside a circle of radius a or
+// with negative roots will be called bad and the rest
+// good. If a is not specified, it will be assumed as 1.
+
+function [goodpoly,badpoly] = polsplit3(fac,a)
+if argn(2) == 1, a = 1; end
+if a>1 error('good polynomial also is unstable'); end
+fac1 = poly(fac(length(fac):-1:1),'z','coeff');
+rts = roots(fac1);
+rts = rts(length(rts):-1:1);
+
+// extract good and bad roots
+badindex = mtlb_find((abs(rts)>=a-1.0e-5)|(real(rts)<-0.05));
+badpoly = coeff(poly(rts(badindex),'z'));
+goodindex = mtlb_find((abs(rts)<a-1.0e-5)&(real(rts)>=-0.05));
+goodpoly = coeff(poly(rts(goodindex),'z'));
+
+// scale by equating the largest terms
+[m,index] = max(abs(fac));
+goodbad = convol(goodpoly,badpoly);
+goodbad = goodbad(length(goodbad):-1:1);
+factor1 = fac(index)/goodbad(index);
+goodpoly = goodpoly * factor1;
+goodpoly = goodpoly(length(goodpoly):-1:1);
+badpoly = badpoly(length(badpoly):-1:1);
+endfunction;
diff --git a/2048/DEPENDENCIES/pp_basic.sci b/2048/DEPENDENCIES/pp_basic.sci new file mode 100644 index 000000000..c8fb0a356 --- /dev/null +++ b/2048/DEPENDENCIES/pp_basic.sci @@ -0,0 +1,26 @@ +// Design of 2-DOF pole placement controller, as discussed in Sec. 9.2.
+// 9.5
+
+// function [Rc,Sc,Tc,gamma] = pp_basic(B,A,k,phi)
+// calculates pole placement controller
+
+
+function [Rc,Sc,Tc,gamm] = pp_basic(B,A,k,phi)
+
+// Setting up and solving Aryabhatta identity
+[Ag,Ab] = polsplit2(A); dAb = length(Ab) - 1;
+[Bg,Bb] = polsplit2(B); dBb = length(Bb) - 1;
+
+[zk,dzk] = zpowk(k);
+
+[N,dN] = polmul(Bb,dBb,zk,dzk);
+dphi = length(phi) - 1;
+
+[S1,dS1,R1,dR1] = xdync(N,dN,Ab,dAb,phi,dphi);
+
+// Determination of control law
+Rc = convol(Bg,R1); Sc = convol(Ag,S1);
+Tc = Ag; gamm = sum(phi)/sum(Bb);
+
+endfunction;
+
diff --git a/2048/DEPENDENCIES/pp_im.sci b/2048/DEPENDENCIES/pp_im.sci new file mode 100644 index 000000000..f1a512fea --- /dev/null +++ b/2048/DEPENDENCIES/pp_im.sci @@ -0,0 +1,25 @@ +// Pole placement controller using internal model principle, as discussed in Sec. 9.4.
+// 9.8
+
+// function [Rc,Sc,Tc,gamma,phit] = pp_im(B,A,k,phi,Delta)
+// Calculates 2-DOF pole placement controller.
+
+function [Rc,Sc,Tc,gamm] = pp_im(B,A,k,phi,Delta)
+
+// Setting up and solving Aryabhatta identity
+[Ag,Ab] = polsplit3(A); dAb = length(Ab) - 1;
+[Bg,Bb] = polsplit3(B); dBb = length(Bb) - 1;
+
+[zk,dzk] = zpowk(k);
+
+[N,dN] = polmul(Bb,dBb,zk,dzk);
+dDelta = length(Delta)-1;
+[D,dD] = polmul(Ab,dAb,Delta,dDelta);
+dphi = length(phi)-1;
+
+[S1,dS1,R1,dR1] = xdync(N,dN,D,dD,phi,dphi);
+
+// Determination of control law
+Rc = convol(Bg,convol(R1,Delta)); Sc = convol(Ag,S1);
+Tc = Ag; gamm = sum(phi)/sum(Bb);
+endfunction;
diff --git a/2048/DEPENDENCIES/pp_im2.sci b/2048/DEPENDENCIES/pp_im2.sci new file mode 100644 index 000000000..61ba58d1b --- /dev/null +++ b/2048/DEPENDENCIES/pp_im2.sci @@ -0,0 +1,31 @@ +// Pole placement controller without intra sample oscillations, as discussed in Sec. 9.5.
+// 9.13
+
+// function [Rc,Sc,Tc,gamma,phit] = pp_im2(B,A,k,phi,Delta,a)
+// 2-DOF PP controller with internal model of Delta and without
+// hidden oscillations
+
+function [Rc,Sc,Tc,gamm,phit] = pp_im2(B,A,k,phi,Delta,a)
+
+if argn(2) == 5, a = 1; end
+dphi = length(phi)-1;
+
+// Setting up and solving Aryabhatta identity
+[Ag,Ab] = polsplit3(A,a); dAb = length(Ab) - 1;
+[Bg,Bb] = polsplit3(B,a); dBb = length(Bb) - 1;
+
+[zk,dzk] = zpowk(k);
+
+[N,dN] = polmul(Bb,dBb,zk,dzk);
+dDelta = length(Delta)-1;
+[D,dD] = polmul(Ab,dAb,Delta,dDelta);
+
+[S1,dS1,R1,dR1] = xdync(N,dN,D,dD,phi,dphi);
+
+// Determination of control law
+Rc = convol(Bg,convol(R1,Delta)); Sc = convol(Ag,S1);
+Tc = Ag; gamm = sum(phi)/sum(Bb);
+
+// Total characteristic polynomial
+phit = convol(phi,convol(Ag,Bg));
+endfunction;
diff --git a/2048/DEPENDENCIES/pp_pid.sci b/2048/DEPENDENCIES/pp_pid.sci new file mode 100644 index 000000000..7071ab671 --- /dev/null +++ b/2048/DEPENDENCIES/pp_pid.sci @@ -0,0 +1,15 @@ +// Solution to Aryabhatta's identity arising in PID controller design, namely Eq. 9.37 on page 363.
+// 9.20
+
+function [Rc,Sc] = pp_pid(B,A,k,phi,Delta)
+
+// Setting up and solving Aryabhatta identity
+dB = length(B) - 1; dA = length(A) - 1;
+[zk,dzk] = zpowk(k);
+[N,dN] = polmul(B,dB,zk,dzk);
+dDelta = length(Delta)-1;
+[D,dD] = polmul(A,dA,Delta,dDelta);
+dphi = length(phi)-1;
+[Sc,dSc,R,dR] = xdync(N,dN,D,dD,phi,dphi);
+Rc = convol(R,Delta);
+endfunction;
diff --git a/2048/DEPENDENCIES/recursion.sci b/2048/DEPENDENCIES/recursion.sci new file mode 100644 index 000000000..ead807bfb --- /dev/null +++ b/2048/DEPENDENCIES/recursion.sci @@ -0,0 +1,35 @@ +// Recursive computation of Ej and Fj
+// 11.1
+
+function [Fj,dFj,Ej,dEj] = recursion(A,dA,C,dC,j)
+Fo = C; dFo = dC;
+Eo = 1; dEo = 0;
+A_z = A(2:dA+1); dA_z = dA-1;
+zi = 1; dzi = 0;
+for i = 1:j-1
+ if (dFo == 0)
+ Fn1 = 0;
+ else
+ Fn1 = Fo(2:(dFo+1));
+ end
+ dFn1 = max(dFo-1,0);
+ Fn2 = -Fo(1)*A_z; dFn2 = dA-1;
+ [Fn,dFn] = poladd(Fn1,dFn1,Fn2,dFn2);
+ zi = convol(zi,[0,1]); dzi = dzi + 1;
+ En2 = Fn(1)*zi; dEn2 = dzi;
+ [En,dEn] = poladd(Eo,dEo,En2,dEn2);
+ Eo = En; Fo = Fn;
+ dEo = dEn; dFo = dFn;
+end
+if (dFo == 0)
+ Fn1 = 0;
+else
+Fn1 = Fo(2:(dFo+1));
+end;
+dFn1 = max(dFo-1,0);
+Fn2 = -Fo(1)*A_z; dFn2 = dA-1;
+[Fn,dFn] = poladd(Fn1,dFn1,Fn2,dFn2);
+Fj = Fn; dFj = dFn;
+Ej = Eo; dEj = dEo;
+endfunction;
+
diff --git a/2048/DEPENDENCIES/spec1.sci b/2048/DEPENDENCIES/spec1.sci new file mode 100644 index 000000000..c749f9695 --- /dev/null +++ b/2048/DEPENDENCIES/spec1.sci @@ -0,0 +1,23 @@ +// Function to implement spectral factorization, as discussed in sec. 13.1.
+// 13.2
+
+function [r,b,rbbr] = spec1(A,dA,B,dB,rho)
+AA = rho * convol(A,flip(A));
+BB = convol(B,flip(B));
+diff1 = dA - dB;
+dBB = 2*dB;
+for i = 1:diff1
+ [BB,dBB] = polmul(BB,dBB,[0 1],1);
+end
+[rbbr,drbbr] = poladd(AA,2*dA,BB,dBB);
+rts = roots(rbbr); // roots in descending order of magnitude
+rts = flip(rts);
+rtsin = rts(dA+1:2*dA);
+b = 1;
+for i = 1:dA,
+ b = convol(b,[1 -rtsin(i)]);
+end
+br = flip(b);
+bbr = convol(b,br);
+r = rbbr(1) / bbr(1);
+endfunction;
diff --git a/2048/DEPENDENCIES/specfac.sci b/2048/DEPENDENCIES/specfac.sci new file mode 100644 index 000000000..1a8ff1332 --- /dev/null +++ b/2048/DEPENDENCIES/specfac.sci @@ -0,0 +1,34 @@ +// Spectral factorization, to solve Eq. 13.47 on page 471.
+// 13.3
+
+// function [r,b,dAFW] = ...
+// specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF)
+// Implements the spectral factorization for use with LQG control
+// design method of Ahlen and Sternard
+
+function [r,b,dAFW] = ...
+ specfac(A,degA,B,degB,rho,V,degV,W,degW,F,degF)
+AFW = convol(A,convol(W,F));
+dAFW = degA + degF + degW;
+AFWWFA = rho * convol(AFW,flip(AFW));
+BV = convol(B,V);
+dBV = degB + degV;
+BVVB = convol(BV,flip(BV));
+diff1 = dAFW - dBV;
+dBVVB = 2*dBV;
+for i = 1:diff1
+ [BVVB,dBVVB] = polmul(BVVB,dBVVB,[0 1],1);
+end
+[rbb,drbb] = poladd(AFWWFA,2*dAFW,BVVB,dBVVB);
+Rbb = polyno(rbb,'z');
+rts = roots(Rbb);
+rtsin = rts(dAFW+1:2*dAFW);
+b = 1;
+for i = 1:dAFW,
+ b = convol(b,[1 -rtsin(i)]);
+end
+b = real(b);
+br = flip(b);
+bbr = convol(b,br);
+r = rbb(1) / bbr(1);
+endfunction;
diff --git a/2048/DEPENDENCIES/tfvar.sci b/2048/DEPENDENCIES/tfvar.sci new file mode 100644 index 000000000..424624d9f --- /dev/null +++ b/2048/DEPENDENCIES/tfvar.sci @@ -0,0 +1,16 @@ +// Cancellation of common factors and determination of covariance
+// 11.7
+
+// function [N,dN,D,dD,yvar] = tfvar(N,dN,D,dD)
+// N and D polynomials in z^{-1} form; discrete case
+
+function [N,dN,D,dD,yvar] = tfvar(N,dN,D,dD)
+
+[N,dN,D,dD] = l2r(N,dN,D,dD);
+N = N/D(1); D = D/D(1);
+LN = length(N); LD = length(D);
+D1 = D;
+if LD<LN, D1 = [D zeros(1,LN-LD)]; dD1 = dD+LN-LD; end
+H = tf(N,D1,1);//TS=1 (sampling time) has been taken constant in tfvar
+yvar = covar_m(H,1);
+endfunction;
diff --git a/2048/DEPENDENCIES/zpowk.sci b/2048/DEPENDENCIES/zpowk.sci new file mode 100644 index 000000000..4b5b43e8d --- /dev/null +++ b/2048/DEPENDENCIES/zpowk.sci @@ -0,0 +1,7 @@ +// Evaluates z^-k.
+// 9.6
+
+function [zk,dzk] = zpowk(k)
+zk = zeros(1,k+1); zk(1,k+1) = 1;
+dzk = k;
+endfunction
|
