links a siti che spiegano scilab

creo un altro blocco-note dedicato a scilab

per i fortunati possessori di Labview >8.0 e' disponibile il gateway Scilab-Labview

regolatore lqr per nxtway

tentativo di trovare il regolatore lqr (wiki) per il nxt di sybille mayer con scilab
attingendo a piene mani (copia-incolla) dal mitico R. Bucher


function [num,den]=tfdata(G)
num=G('num');
den=G('den');
endfunction

A=[0,1,0,0; 94.05,0,0,1.3;0,0,0,1; -544.05,0,0,-14.3]
B=[0;-0.21;0;2.32]
C=[1,0,0,0;0,0,1,0]
D=[0;0]
nxtway_SYB=syslin('c',A,B,C,D)
s=%s;
Ts=0.03
spec(A)
autovalori_di_A=spec(A)
sys=syslin('c',A,B,C,D);
sysd=dscr(sys,Ts);
[ad,bd,cd,dd]=abcd(sysd);

// pesi del controllore LQR

Q=diag([10000,1,10000,1]); // 4 per 4
R=[1]; // 1 per 1

//calcola i guadagni LQR per un sistema discreto
[n1,d1]=size(ad);
big=sysdiag(Q,R);
[w,wp]=fullrf(big ,1e-20);
C1=wp(:,1:n1);
D12=wp(:,n1+1:$);
P=syslin('d',ad,bd,C1,D12);
[k,X]=lqr(P)
k_lqr= -k;

E=spec(ad);
plqr=abs(real(log(E)/Ts))';
pmax=max(plqr);
k_lqr=-k_lqr;
preg=spec(A);
// osservatore di ordine ridotto
poli_oss=exp([real(preg(3)),real(preg(4))]*10*Ts);
T=[0,0,0,1;0,1,0,0];

// trova l'ossrvatore di ordine ridotto per A,B,C,D, con polo di osservatore
// T e' la matrice necessaria per rendere [C;T] invertibile
P=[cd;T]
invP=inv([cd;T])
AA=P * ad * invP
ny=size(cd,1)
nx=size(ad,1)
nu=size(bd,2)
A11=AA(1:ny,1:ny)
A12=AA(1:ny,ny+1:nx)
A21=AA(ny+1:nx,1:ny)
A22=AA(ny+1:nx,ny+1:nx)
L1=ppol (A22',A12',poli_oss)';
nn=nx - ny;
A_redobs=[-L1 eye(nn,nn)]*P*ad*invP*[zeros(ny,nn); eye(nn,nn)];
B_redobs=[-L1 eye(nn,nn)]*[P*bd P*ad*invP*[eye(ny,ny);L1]]*[eye(nu,nu) zeros(nu,ny);-dd, eye(ny,ny)];
C_redobs=invP*[zeros(ny,nx-ny);eye(nn,nn)];
D_redobs=invP*[zeros(ny,nu) eye(ny,ny);zeros(nx-ny,nu) L1]*[eye(nu,nu) zeros(nu,ny);-dd, eye(ny,ny)];
ao=A_redobs
bo=B_redobs
co=C_redobs
do=D_redobs

// Crea la forma compatta dell'osservatore ABCD e il guadagno K
//
//

Bu=bo(:,1);
By=bo(:,2:$);
Du=do(:,1);
Dy=do(:,2:$);
X=inv(1+k_lqr*Du);
Ac=ao - Bu*X*k_lqr*co;
Bc=[Bu*X,By-Bu*X*k_lqr*Dy]
Cc=-X*k_lqr*co;
Dc=[X,-X*k_lqr*Dy]
Greg=syslin('d',Ac,Bc,Cc,Dc)
///////////////////////////////////////////
Gregtf=ss2tf(Greg);
Gregtf( :,2)
Gregtf( :,3)
[g1n,g1d]=tfdata(Gregtf( :,2))
[g2n,g2d]=tfdata(Gregtf( :,3))


questo butta fuori come regolatore il seguente:

Greg =


Greg(1) (state-space system:)

!lss A B C D X0 dt !

Greg(2) = A matrix =

- 1.1704196 - 16.638285
0.1322204 2.048937

Greg(3) = B matrix =

0.0329866 - 559.33453 - 69.773154
- 0.0036029 41.850078 4.3785807

Greg(4) = C matrix =

- 35.525535 - 497.8127

Greg(5) = D matrix =

1. - 16487.356 - 1308.8181

Greg(6) = X0 (initial state) =

0.
0.

Greg(7) = Time domain =

d

regolatore con scilab

Riferimento al problema di regolare il sistema instabile del nxt-way di Sybille Mayer.

Con scilab:

A=[0,1,0,0; 94.05,0,0,1.3;0,0,0,1; -544.05,0,0,-14.3] // pag 24
B=[0;-0.21;0;2.32]
C=[1,0,0,0;0,0,1,0]
D=[0;0]
nxtway_SYB=syslin('c',A,B,C,D) // genera sistema lineare
s=%s; // variabile di laplace
spec(A) // calcola autovalori
autovalori_di_A=spec(A)
poli_desiderati=[-2.6 -3.9 -5.2 -6.5] // poli desiderati
K=ppol(A,B,poli_desiderati) // calcolo della matrice di feedback K

F=A-B*K // sistema in retroazione
g2=syslin('c',F,B,C,D) // genera nuovo sistema lineare
spec(F)
autovalori_di_F=spec(F) //calcola autovalori, per verificare il piazzamento poli

si ottengono i seguenti risultati:

--> A =

0. 1. 0. 0.
94.05 0. 0. 1.3
0. 0. 0. 1.
- 544.05 0. 0. - 14.3
B =

0.
- 0.21
0.
2.32
C =

1. 0. 0. 0.
0. 0. 1. 0.
D =

0.
0.
nxtway_SYB =


nxtway_SYB(1) (state-space system:)

!lss A B C D X0 dt !

nxtway_SYB(2) = A matrix =

0. 1. 0. 0.
94.05 0. 0. 1.3
0. 0. 0. 1.
- 544.05 0. 0. - 14.3

nxtway_SYB(3) = B matrix =

0.
- 0.21
0.
2.32

nxtway_SYB(4) = C matrix =

1. 0. 0. 0.
0. 0. 1. 0.

nxtway_SYB(5) = D matrix =

0.
0.

nxtway_SYB(6) = X0 (initial state) =

0.
0.
0.
0.

nxtway_SYB(7) = Time domain =

c
ans =

0
7.8847523
- 4.598569
- 17.586183
autovalori_di_A =

0
7.8847523
- 4.598569
- 17.586183
poli_desiderati =

- 2.6 - 3.9 - 5.2 - 6.5
K =

- 1063.3266 - 123.77134 - 3.2972279 - 9.5224059
F =

0. 1. 0. 0.
- 129.2486 - 25.991982 - 0.6924179 - 0.6997052
0. 0. 0. 1.
1922.8678 287.14951 7.6495687 7.7919818
g2 =


g2(1) (state-space system:)

!lss A B C D X0 dt !

g2(2) = A matrix =

0. 1. 0. 0.
- 129.2486 - 25.991982 - 0.6924179 - 0.6997052
0. 0. 0. 1.
1922.8678 287.14951 7.6495687 7.7919818

g2(3) = B matrix =

0.
- 0.21
0.
2.32

g2(4) = C matrix =

1. 0. 0. 0.
0. 0. 1. 0.

g2(5) = D matrix =

0.
0.

g2(6) = X0 (initial state) =

0.
0.
0.
0.

g2(7) = Time domain =

c
ans =

- 6.5
- 2.6
- 3.9
- 5.2
autovalori_di_F =

- 6.5
- 2.6
- 3.9
- 5.2


si vede che il sistema retroazionato ha i poli (negativi) del valore che si voleva.

I poli ottenuti da Sybille sono:
[-1051.93 -122.34 -3.32 -9.44]
mentre lo script scilab da'
K = - 1063.3266 - 123.77134 - 3.2972279 - 9.5224059

piu' o meno gli stessi.

risposta all'impulso instabile, nxtway_SYB






















risposta all'impulso stabile, g2























risposta al gradino sistema instabile, nxtway_SYB:

























risposta al gradino sistema stabile, g2:

lego NXT segway

(immagine png, scaricabile e zoomabile)

questo e' il programma NXT-G che realizza il PID (numeri interi) per il NXT-way standard, che in nbc viene implementato cosi':

http://www.philohome.com/nxtway/nxtway.nbc

(da http://www.philohome.com/nxtway/nxtway.htm)

il programma touchway5.rbt in figura e' stato realizzato da

http://www.freewebs.com/thesystemprogrammer/


per approfondimenti, vedere tesi di Serra

http://www.epokh.org/tesi/TesiMathSerra.pdf

E

http://nikemagic.altervista.org/download/Servo/progetto_pendolo_pid.pdf

http://www.plcforum.info/didattica/conreg/conreg.htm

E, ma in inglese,

http://www.microchip.com/stellent/idcplg?IdcService=SS_GET_PAGE&nodeId=1824&appnote=en021807















































da http://www.pages.drexel.edu/~ttl28/thesis.html#labview :
Modeling/Simulation

One major challange is how to properly model this system i.e. a specific transfer function. To identify the system, I first began but re-designing the model in SolidWorks to find the robot's Center-of-Mass and Moment-of-Inertia to use in plant designs of mobile inverted pendulums.

qui simulazioni con matlab

invece, per i poveracci, le simulazioni si possono fare con scilab, seguendo l'esempio:

utile guida per scicos/scilab : Bucher

eserciziario

sistemi MIMO in scilab: http://www.wolffdata.se/scilab/mimo1.html

varia didattica

di interesse:

http://madrobotics.blogspot.com/

Nel caso si abbiano 2 sensori di luce a disposizione (in realta' nel testo si dice che funziona anche con uno, quello dietro, "hinter"):

http://www.uni-kassel.de/fb16/rat/diplom/Diplom1Mayer.pdf

tale tesi ad opera di Sybille Mayer, in tedesco, e' veramente interessante.

Viene presentato il modello di pendolo inverso IMPERNIATO SU RUOTE (cosa non affrontata in altri lavori); dopo linearizzazione viene identificato il modello del motore
(risposta al gradino) e il modello del sensore di luce (con interessante digressione sul comportamento su pavimenti di colore diverso...).

Il modello meccanico e' a pagina 24.

A questo punto SM progetta il regolatore sfruttando sia il sensore di luce (anche 1 solo) e il sensore di rotazione nel motore. L'idea e' un PD sulla luce,cioe' l'angolo alfa, e un PID sulla rotazione cioe' l'angolo gamma. Per intenderci alfa+gamma e' la rotazione della ruota, che moltiplicata per il raggio da' lo spostamento del NXTway dalla posizione iniziale. Il sistema totale e' a pagina 36 (PD+PD), con variante (PD+PID) alla pagina 37.

Le simulazioni sono basate su questo modello del sistema (simulink):

Il programma per NXT-G (sic!) e' scaricabile qui:

http://forums.nxtasy.org/index.php?act=attach&type=post&id=616

NB: nel programma i coefficienti sono leggermente diversi da quelli teorici.

codice scilab per analizzare il sistema di Sybille da lanciare con exec nomefile:

A=[0,1,0,0; 94.05,0,0,1.3;0,0,0,1; -544.05,0,0,-14.3];
B=[0;-0.21;0;2.32];
C=[1,0,0,0;0,0,1,0];
D=[0;0];
nxtway_SYB=syslin('c',A,B,C,D);
s=%s;
nxtway_SYBtf=ss2tf(nxtway_SYB) // passo nel dominio Laplace
cl_nxtway_SYBtf=clean(nxtway_SYBtf) // mette a zero i coefficienti infinitesimi
t=0:0.01:1;
u=ones(1:length(t)); // step input
y=csim(u,t,nxtway_SYB); // risposta al gradino instabile
xset('window',1);
xbasc()
plot2d(t,y(2,:))
plot2d(t,y(1,:), style=5)

// poli del sistema
spec(A)
autovalori_di_A=spec(A)

// dalle specifiche ai poli desiderati
os=9.48 // overshoot desiderato
ts=0.74 // settling time desiderato
os=os/100
xi=-log(os)/sqrt(%pi*%pi+log(os)*log(os)) //smorzamento xi
tetaxi=acos(xi)
wn=(-log(0.02*sqrt(1-xi*xi)))/(ts*xi) // risonanza wn
th=acos(xi)
ps=-xi*wn+%i*wn*sqrt(1-xi*xi); // uno dei 2 poli complessi coniugati del sys 2o ordine per fare os e ts
cl_nxtway_SYBtf(2,1)
trfmod(cl_nxtway_SYBtf(2,1),'p') // calcolo dei poli e degli zeri
poli_desiderati=[-6 -2 ps conj(ps)] // poli sono i desiderati
K=ppol(A,B,poli_desiderati) // calcolo della matrice di feedback K
F=A-B*K // sistema in retroazione
g2=syslin('c',F,B,C,D)
t=0:0.01:3;
u=ones(1:length(t)); // step input
y=csim(u,t,g2); // risposta al gradino controllato
xset('window',1);
xbasc()
plot2d(t,y(2,:))
plot2d(t,y(1,:), style=5)


scicos

Altro progetto, in tedesco, con NXT controllato voa USB (purtroppo Labview8.2) oppure standalone:

http://projekte6.fhnw.ch/technik/eit/Herbst2007/BruWid/

In italiano:

scienzaludica.it

Infine una bella carrellata di robot affini:
http://leiwww.epfl.ch/joe/joerelatives.html

interessante pendolo pilotato verticalmente:

http://www.zfm.ethz.ch/~leine/vertically_driven_pendulum.htm

con istruzioni

http://www.zfm.ethz.ch/~leine/LEGO/VerticallyDrivenPendulum/VDPBuildinginstructions.pdf

anche questo:

http://www.zehl.com/?nxtway