[PYTHON] Introduction d'une méthode de conception de contrôleur pilotée par les données

introduction

Le contrôle de rétroaction a été utilisé pour «contrôler» tout, des appareils ménagers aux équipements industriels. Ceux qui effectuent ce contrôle sont appelés contrôleurs et diverses méthodes de conception ont été envisagées. Cette fois, je voudrais présenter une approche relativement nouvelle appelée «méthode de conception de contrôleur pilotée par les données».

De plus, je laisserai la partie théorique au papier et écrirai principalement sur son utilisation. Donc, je voudrais l'expliquer avec le code de MATLAB et Python.

Du contrôle piloté par les données

En gros, il existe trois approches de la conception des contrôleurs:

1. Réglage manuel des paramètres (méthode de sensibilité limite, etc.)
2. Conception de contrôleur basée sur un modèle (régulateur optimal, etc.)
3. Autres ($ H_ {\ infty} $ control, contrôle piloté par les données, etc.)

Puisque 1 est manuel, il n'est pas difficile d'imaginer que la conception est difficile car elle nécessite des essais et des erreurs. Par conséquent, comme en 2, une méthode a été proposée dans laquelle le * modèle mathématique à contrôler a été exactement obtenu * et le modèle a été utilisé pour la conception. La tâche de trouver exactement ce modèle est appelée «identification». En décrivant le problème d'optimisation à l'aide du modèle identifié, nous sommes maintenant en mesure d'obtenir le contrôleur théoriquement optimal. Cependant, c'est le musicien qui veut exactement ce ** modèle **, et il faut beaucoup de temps et d'efforts pour concilier théorie et expérience. Par conséquent, d'autres méthodes fondées sur l'hypothèse que ** les modèles ne peuvent pas être obtenus exactement ** ont été envisagées. Le contrôle basé sur les données en fait partie.

Qu'est-ce que le contrôle piloté par les données?

Selon le site Web du professeur Kaneko de l'Université des télécommunications, qu'est-ce que le contrôle piloté par les données?

Une méthode de conception d'un contrôleur directement à partir des données et du comportement de l'appareil à contrôler lorsqu'il est piloté.

est. En d'autres termes, le contrôleur peut être conçu sans créer de modèle, tant qu'il y a des données de signal pendant la conduite.

Il existe actuellement deux méthodes de conception de contrôleurs basées sur les données bien connues [^ excuse IFT].

[^ excuseIFT]: En fait, historiquement, le contrôle non-preuve et l'IFT peuvent avoir besoin d'explications, mais je les omettrai.

Cette fois, j'ai écrit le code sur VRFT afin que le problème d'optimisation puisse être résolu facilement.

Exemple d'environnement d'exécution de code

• MATLAB 2017b
  • Control System Toolbox
  • System Identification Toolbox
• Python 3.7.5 (environnement Anaconda)
  • matplotlib
  • numpy
  • scipy
  • python-control

Notez que le code Python dans l'explication omet l'instruction d'importation. Voir [Code entier](#Whole Code) pour le code de travail.

Procédure VRFT

Il n'y a que trois choses qu'un concepteur doit faire pour utiliser VRFT.

1. Détermination des spécifications de conception
2. Acquisition des données d'entrée / sortie
3. Résoudre le problème d'optimisation selon l'algorithme

Cette fois, je vais expliquer la conception du contrôleur du système de contrôle de vitesse en utilisant un moteur comme exemple.

Détermination des spécifications de conception

Tout d'abord, décidons des spécifications de conception. Les trois spécifications de conception suivantes peuvent être déterminées.

Par exemple, si vous voulez que le contrôleur PID fournisse une réponse avec une constante de temps $ \ tau = 0,02 $ [s],

\begin{align}
M(s) &= \frac{1}{0.02 s + 1} \\&= \frac{50}{s + 50}\\
\boldsymbol{\beta}(s) &= \begin{bmatrix}1& \frac{1}{s}& s\end{bmatrix}^\mathrm{T}
\end{align}

Réglez comme suit. De plus, étant donné qu'il existe une forte possibilité que le bruit soit mélangé dans les données dans la gamme de hautes fréquences, il suffit souvent de régler un filtre passe-bas pour le poids de fréquence.

Le code ci-dessus est le suivant.

vrft.m

%Temps d'échantillonnage et opérateurs
Ts = 0.001;
s = tf('s');        %Opérateurs de Laplace
z = tf('z', Ts);    %Opérateur d'avance temporelle z

%Modèle de référence M
tau = 0.02;
M = c2d(1/(tau*s + 1), Ts); %Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

%Fonction poids
gW = 100;
W = c2d(gW/(s + gW), Ts);   %Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

%Structure de contrôle beta
beta = minreal([1; Ts/(1 - z^-1); (1 - z^-1)/Ts]); %Contrôleur PID

vrft.py

# ===Détermination des spécifications de conception===
#Temps d'échantillonnage
Ts = 0.001

#Modèle de référence M
tau = 0.02
M = ctl.c2d(ctl.tf([1], [tau, 1]), Ts)  #Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

#Fonction poids W
gW = 100
W = ctl.c2d(ctl.tf([gW], [1, gW]), Ts)  #Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

#Structure de contrôle beta
A, B, C, D = abcd_normalize(
    [[1., 0], [0, 0]],
    [[1.], [-1.]],
    [[0, 0], [Ts, 0], [0, -1/Ts]],
    [[1.], [0], [1/Ts]])
beta = ctl.ss(A, B, C, D, Ts)
del A, B, C, D

Acquisition des données d'entrée / sortie

Ensuite, les données d'entrée / sortie $ \ {u_0 (t), y_0 (t) } $ sont acquises à partir de la cible de contrôle. Dans le cas de VRFT, les données sont acquises afin que les conditions suivantes soient remplies.

• Assurez-vous d'obtenir la réponse de la boucle ouverte (seulement $ P (s) $ sans contrôle)
• Si possible, utilisez le signal d'entrée blanc $ u_0 (t) $

En tant que signal de blancheur, il est recommandé d'utiliser le signal de la série M qui est souvent utilisé pour l'identification du système. Dans le cas de MATLAB, les signaux de la série M peuvent être facilement générés à l'aide de la fonction idinput de la boîte à outils d'identification du système. Dans le cas de Python, il peut être généré par la [fonction max_len_seq] de Scipy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.max_len_seq.html?highlight=correlation). Cependant, il semble que des signaux d'une série différente de MATLAB soient générés, donc dans cette version de Python, j'ai créé ma propre fonction prbs () et je l'ai utilisée.

Ici, un modèle de l'objet contrôlé (moteur dans ce cas) est nécessaire pour effectuer la simulation [^ note1]. Donc, à l'origine, j'ai pensé qu'il était nécessaire d'expliquer le modèle du moteur, mais Manao's Great article J'ai trouvé des éléments / ed383a30bbe0c7b9534a). Par conséquent, nous allons omettre la modélisation ici et supposer que le moteur peut être représenté par un système primaire.

[^ note1]: Notez que cela n'est nécessaire que pour la récupération de données et qu'aucun modèle n'est utilisé dans la conception. (En fait, il peut être conçu comme «clearvars P» après l'acquisition de données.)

Par conséquent, dans la simulation, les données peuvent être acquises avec le code suivant.

vrft.m

%Signal d'entrée utilisé pour l'acquisition de données u(signal série m)
n = 15;                                     %Nombre d'étapes
T = 2^n - 1;                                %Nombre de données par cycle
p = 15;                                     %Nombre de cycles de données
N = T*p;                                    %Le nombre de données
u0 = idinput([T 1 p],'prbs',[0,1],[-1,1]);   %Vecteur de signal

%Densité spectrale de puissance du signal d'entrée phi_u
phi_u = 1;                  %Le signal d'entrée est supposé être blanc

%Modèle contrôlé P
Tp = 0.74;
Kp = 1.02;
P = c2d(Kp/(Tp*s + 1), Ts);

%Signal de sortie y0
y0 = lsim(P, u0);

vrft.py

# ===Acquisition des données d'entrée / sortie===
def prbs(n, p):
    # matlab compatible PRBS signal generator

    taps = {3: [0, -1], 4: [0, -1], 5: [1, -1], 6: [0, -1], 7: [0, -1],
            8: [0, 1, 6, -1], 9: [3, -1], 10: [2, -1], 11: [8, -1],
            12: [5, 7, 10, -1], 13: [8, 9, 11, -1], 14: [3, 7, 12, -1],
            15: [13, -1], 16: [3, 12, 14, -1], 17: [13, -1], 18: [10, -1]}
    N = (2**n - 1)*p
    x = np.ones(n, dtype=np.int8)
    u = np.zeros(N, dtype=np.int8)
    tap = taps[n]
    for i in range(N):
        u[i] = x[-1]
        x0 = x[tap[0]] ^ x[tap[1]]
        x = np.roll(x, 1)
        x[0] = x0
    return u


#Signal d'entrée utilisé pour l'acquisition de données u(signal série m)
n = 15
T = 2**n - 1
p = 15
N = T*p
# u0, _ = max_len_seq(n, length=N)  #La série de signaux est différente de celle de MATLAB
u0 = prbs(n, p)
u0 = -(2.*u0 - 1.)

#Densité spectrale de puissance du signal d'entrée phi_u
phi_u = 1.

#Modèle contrôlé P
Tp = 0.74
Kp = 1.02
P = ctl.c2d(ctl.tf([Kp], [Tp, 1]), Ts)

#Signal de sortie y0
y0, t0, _ = ctl.lsim(P, u0)

Résoudre le calcul d'optimisation selon l'algorithme

Enfin, le signal est généré selon l'algorithme VRFT et le calcul d'optimisation est résolu.

C'est un problème de minimisation pour la fonction d'évaluation $ J_ \ mathrm {MR} (\ boldsymbol {\ rho}) $, car VRFT vise à répondre au plus près du modèle de référence $ M (z) $.

\begin{align}
J_\mathrm{MR}(\boldsymbol{\rho}) &= \left\| \left( \frac{P(z) C(z, \boldsymbol{\rho})}{1 + P(z) C(z, \boldsymbol{\rho})} - M(z) \right) W(z) \right\|_2 ^2
\end{align}

Ici, $ \ boldsymbol {\ rho} $ représente le paramètre de contrôle, $ C (z, \ boldsymbol {\ rho}) = \ boldsymbol {\ rho} ^ {\ rm T} \ boldsymbol {\ beta } (z) Il existe une relation $. Cependant, puisque cette expression utilise le modèle $ P (z) $ pour être contrôlée, il est nécessaire de la réécrire avec des données d'entrée / sortie. L'explication théorique est coupée, mais le problème de la minimisation de la fonction d'évaluation $ J_ \ mathrm {MR} $ est le suivant, et le problème de la minimisation de $ J_ \ mathrm {VR} ^ N $ est graduel. Est équivalent à.

\begin{align}
J_{\rm VR}^{N}(\boldsymbol{\rho}) &= \frac{1}{N} \sum_{t=1}^{N}(u_L(k) - \boldsymbol{\rho}^{\rm T} \boldsymbol{\varphi}(t))^2
\end{align}

Ici, les données de séries temporelles $ u_L (t) $ et $ \ boldsymbol {\ rho} $ peuvent être générées par la formule suivante en utilisant le filtre $ L (z) $.

\begin{align}
u_L(t) &= L(z) u_0(t)\\
y_L(t) &= L(z) y_0(t)\\
%\boldsymbol{\varphi}(k) &= \boldsymbol{\beta}(z) (\tilde{r}_L(t) - y_L(t))
e_L(t) &= (M^{-1} - 1)y_L(t)\\
\boldsymbol{\varphi}(k) &= \boldsymbol{\beta}(z) e_L(t)
\end{align}

Cependant, le filtre $ L (z) $ doit satisfaire ce qui suit.

\begin{align}
|L(z)|^2 &= \frac{|1-M(z)|^2 |M(z)|^2 |W(z)|^2}{\phi_u(\omega)} &\forall \omega
\end{align}

Ensuite, $ J_ \ mathrm {VR} ^ N (\ boldsymbol {\ rho}) $ n'utilise pas le modèle contrôlé $ P (z) $, et les données d'entrée / sortie $ \ {u_0 (t), y_0 (t) Il ne peut être décrit qu'avec) } $. De plus, puisque cette équation est linéaire par rapport au paramètre de contrôle $ \ boldsymbol {\ rho} $, elle peut être résolue par la méthode des moindres carrés.

D'après ce qui précède, la partie qui résout le problème d'optimisation est le code suivant.

vrft.m

%Préfiltre L
L = minreal(M*(1 - M)/phi_u);

%Signal d'entrée filtré ul
ul = lsim(L, u0);

%Signal de pseudo-erreur el
el = lsim(L*(M^(-1) - 1), y0);

%Contrôle de la sortie avant les paramètres phi
phi = lsim(beta, el);

%Paramètre optimal rho
rho = phi\ul;               %Résolvez la méthode des moindres carrés sous forme de matrice(mldivide)

%Contrôleur conçu C
C = minreal(rho.' * beta);  %Demandez un contrôleur

%Fonction d'évaluation Jmr
Jmr = mean(ul - phi * rho); %Vérifiez la fonction d'évaluation au format matriciel

vrft.py

# ===Conçu par VRFT===
#Préfiltre L
L = ctl.minreal(M*(1 - M)/phi_u)

#Signal d'entrée filtré ul
ul, _, _ = ctl.lsim(L, u0)

#Signal de pseudo-erreur el
el, _, _ = ctl.lsim(ctl.minreal(L*(M**-1 - 1)), y0.flatten())

#Contrôle de la sortie avant les paramètres
phi, _, _ = ctl.lsim(beta, el.flatten())

#Paramètre optimal rho
solution = np.linalg.lstsq(phi, ul, rcond=None)
rho = solution[0]

#Contrôleur conçu C
C = ctl.ss([0], [0, 0, 0], [0], rho.T, Ts) * beta   #Demandez un contrôleur

#Fonction d'évaluation Jmr
Jmr = np.mean(ul - np.dot(phi, rho))    #Vérifiez la fonction d'évaluation au format matriciel

Confirmation des performances

Avec la procédure ci-dessus, il semble que le contrôleur pourrait être conçu en utilisant uniquement les données d'entrée / sortie sans utiliser le modèle à contrôler. Les paramètres de contrôle optimaux $ \ rho $ obtenus sont les suivants. Vous pouvez voir que MATLAB et Python donnent presque les mêmes résultats.

[^ ideal]: $ M (z) = \ frac {C ^ * (z) P (z)} {1 + C ^ * (z) P (z)} Un contrôleur qui satisfait $ C ^ * $ Valeur (c'est ce qu'on appelle le contrôleur idéal)

Maintenant, vérifions si la réponse du modèle de référence est réellement réalisée par ce contrôle. Le code à vérifier est le suivant.

vrft.m

%Système complet avec contrôleur G
G = minreal(feedback(P*C, 1));

%Réponse pas à pas
fig1 = figure('name', 'Step plot');
stepplot(G, M);
legend({'$$\frac{CP}{1+CP}$$', '$M$'},...
    'Interpreter', 'latex', 'FontSize', 14, 'Location', 'southeast');

%Affichage du diagramme de la carte
fig2 = figure('name', 'Bode plot of controller');
bodeplot(G, M, {1,100});

vrft.py

# ===Confirmation des performances===
#Système complet avec contrôleur G
G = ctl.minreal(ctl.feedback(P*C, 1))

#Réponse pas à pas
plt.figure()
plt.title("Step response of closed loop system")
timeRange = np.arange(0, 0.12 + Ts, Ts)
ym, _ = ctl.step(M, timeRange)
yg, _ = ctl.step(G, timeRange)
plt.plot(timeRange, ym, timeRange, yg)
plt.xlabel("Time [s]", usetex=True)
plt.ylabel("Velocity [V]")
plt.legend(['Reference model', 'Closed loop system'], loc='lower right')

#Affichage du diagramme de la carte
plt.figure()
plt.title("Bode plot of closed loop system")
ctl.bode(M, G)
plt.legend(['Reference model', 'Closed loop system'], loc='lower left')
plt.show()

Ici, seule la réponse échelonnée est affichée sous forme d'image. Vous pouvez voir que la réponse du modèle de référence peut être reproduite avec une précision indiscernable.

Code entier

Enfin, je vais vous montrer tout le code que j'ai introduit.

%Script de conception VRFT
%Conception pour le système de contrôle de vitesse du moteur.

% Copyright (c) 2019 larking95(https://qiita.com/larking95)
% Released under the MIT Licence 
% https://opensource.org/licenses/mit-license.php

%%Initialisation
clearvars;
close all;

%%Détermination des spécifications de conception
%Temps d'échantillonnage et opérateurs
Ts = 0.001;
s = tf('s');        %Opérateurs de Laplace
z = tf('z', Ts);    %Opérateur d'avance temporelle z

%Modèle de référence M
tau = 0.02;
M = c2d(1/(tau*s + 1), Ts); %Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

%Fonction poids
gW = 100;
W = c2d(gW/(s + gW), Ts);   %Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

%Structure de contrôle beta
beta = minreal([1; Ts/(1 - z^-1); (1 - z^-1)/Ts]); %Contrôleur PID

%%Acquisition des données d'entrée / sortie
%Signal d'entrée utilisé pour l'acquisition de données u(signal série m)
n = 15;                                     %Nombre d'étapes
T = 2^n - 1;                                %Nombre de données par cycle
p = 15;                                     %Nombre de cycles de données
N = T*p;                                    %Le nombre de données
u0 = idinput([T 1 p],'prbs',[0,1],[-1,1]);   %Vecteur de signal

%Densité spectrale de puissance du signal d'entrée phi_u
phi_u = 1;                  %Le signal d'entrée est supposé être blanc

%Modèle contrôlé P
Tp = 0.74;
Kp = 1.02;
P = c2d(Kp/(Tp*s + 1), Ts);

%Signal de sortie y0
y0 = lsim(P, u0);

%%Conçu par VRFT
%Préfiltre L
L = minreal(M*(1 - M)/phi_u);

%Signal d'entrée filtré ul
ul = lsim(L, u0);

%Signal de pseudo-erreur el= 
el = lsim(L*(M^(-1) - 1), y0);

%Contrôle de la sortie avant les paramètres phi
phi = lsim(beta, el);

%Paramètre optimal rho
rho = phi\ul;               %Résolvez la méthode des moindres carrés sous forme de matrice(mldivide)

%Contrôleur conçu C
C = minreal(rho.' * beta);  %Demandez un contrôleur

%Fonction d'évaluation Jmr
Jmr = mean(ul - phi * rho); %Vérifiez la fonction d'évaluation au format matriciel

%%Confirmation des performances
%Système complet avec contrôleur G
G = minreal(feedback(P*C, 1));

%Réponse pas à pas
fig1 = figure('name', 'Step plot');
stepplot(G, M);

%Affichage du diagramme de la carte
fig2 = figure('name', 'Bode plot of controller');
bodeplot(G, M, {1,100});

# -*- coding: utf-8 -*-
"""
Script de conception VRFT
Conception pour le système de contrôle de vitesse du moteur.

Copyright (c) 2019 larking95(https://qiita.com/larking95)
Released under the MIT Licence 
https://opensource.org/licenses/mit-license.php
"""

import matplotlib.pyplot as plt
import numpy as np
import control.matlab as ctl

# from scipy.signal import max_len_seq
from scipy.signal import abcd_normalize


def prbs(n, p):
    # matlab compatible PRBS signal generator

    taps = {3: [0, -1], 4: [0, -1], 5: [1, -1], 6: [0, -1], 7: [0, -1],
            8: [0, 1, 6, -1], 9: [3, -1], 10: [2, -1], 11: [8, -1],
            12: [5, 7, 10, -1], 13: [8, 9, 11, -1], 14: [3, 7, 12, -1],
            15: [13, -1], 16: [3, 12, 14, -1], 17: [13, -1], 18: [10, -1]}
    N = (2**n - 1)*p
    x = np.ones(n, dtype=np.int8)
    u = np.zeros(N, dtype=np.int8)
    tap = taps[n]
    for i in range(N):
        u[i] = x[-1]
        x0 = x[tap[0]] ^ x[tap[1]]
        x = np.roll(x, 1)
        x[0] = x0
    return u


# ===Détermination des spécifications de conception===
#Temps d'échantillonnage
Ts = 0.001

#Modèle de référence M
tau = 0.02
M = ctl.c2d(ctl.tf([1], [tau, 1]), Ts)  #Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

#Fonction poids W
gW = 100
W = ctl.c2d(ctl.tf([gW], [1, gW]), Ts)  #Système de décalage du premier ordre discrétisé par maintien d'ordre zéro

#Structure de contrôle beta
A, B, C, D = abcd_normalize(
    [[1., 0], [0, 0]],
    [[1.], [-1.]],
    [[0, 0], [Ts, 0], [0, -1/Ts]],
    [[1.], [0], [1/Ts]])
beta = ctl.ss(A, B, C, D, Ts)
del A, B, C, D

# ===Acquisition des données d'entrée / sortie===
#Signal d'entrée utilisé pour l'acquisition de données u(signal série m)
n = 15
T = 2**n - 1
p = 15
N = T*p
# u0, _ = max_len_seq(n, length=N)  #La série de signaux est différente de celle de MATLAB
u0 = prbs(n, p)
u0 = -(2.*u0 - 1.)

#Densité spectrale de puissance du signal d'entrée phi_u
phi_u = 1.

#Modèle contrôlé P
Tp = 0.74
Kp = 1.02
P = ctl.c2d(ctl.tf([Kp], [Tp, 1]), Ts)

#Signal de sortie y0
y0, t0, _ = ctl.lsim(P, u0)

# ===Conçu par VRFT===
#Préfiltre L
L = ctl.minreal(M*(1 - M)/phi_u)

#Signal d'entrée filtré ul
ul, _, _ = ctl.lsim(L, u0)

#Signal de pseudo-erreur el
el, _, _ = ctl.lsim(ctl.minreal(L*(M**-1 - 1)), y0.flatten())

#Contrôle de la sortie avant les paramètres
phi, _, _ = ctl.lsim(beta, el.flatten())

#Paramètre optimal rho
solution = np.linalg.lstsq(phi, ul, rcond=None)
rho = solution[0]

#Contrôleur conçu C
C = ctl.ss([0], [0, 0, 0], [0], rho.T, Ts) * beta   #Demandez un contrôleur

#Fonction d'évaluation Jmr
Jmr = np.mean(ul - np.dot(phi, rho))    #Vérifiez la fonction d'évaluation au format matriciel

# ===Confirmation des performances===
#Système complet avec contrôleur G
G = ctl.minreal(ctl.feedback(P*C, 1))

#Réponse pas à pas
plt.figure()
plt.title("Step response of closed loop system")
timeRange = np.arange(0, 0.12 + Ts, Ts)
ym, _ = ctl.step(M, timeRange)
yg, _ = ctl.step(G, timeRange)
plt.plot(timeRange, ym, timeRange, yg)
plt.xlabel("Time [s]", usetex=True)
plt.ylabel("Velocity [V]")
plt.legend(['Reference model', 'Closed loop system'], loc='lower right')

#Affichage du diagramme de la carte
plt.figure()
plt.title("Bode plot of closed loop system")
ctl.bode(M, G)
plt.legend(['Reference model', 'Closed loop system'], loc='lower left')
plt.show()

en conclusion

Cette fois, tout en omettant complètement la théorie, j'ai présenté les grandes lignes de la conception du contrôleur piloté par les données et le code réel de l'un d'entre eux, VRFT. Ingénierie du contrôle Je n'en connais pas autant que les autres participants au calendrier de l'Avent, mais j'espère que cela aidera tout le monde à s'intéresser au contrôle basé sur les données.

De plus, si vous avez des inquiétudes, n'hésitez pas à commenter.

Les références

1. Introduction de VRFT (site de Campi et Savaresi)
2. Introduction of FRIT (article d'introduction par le professeur Kaneko)
3. [Principes de base de l'identification du système](https://www.amazon.co.jp/%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0%E5%90% 8C% E5% AE% 9A% E3% 81% AE% E5% 9F% BA% E7% A4% 8E-% E8% B6% B3% E7% AB% 8B-% E4% BF% AE% E4% B8% 80 / dp / 4501114800) (Livre du Dr Adachi)
4. Identification système du moteur à courant continu utilisant Arduino et MATLAB / Simulink (Manao article )
5. NumPy for Matlab users(Scipy.org)