7.4.4 Extender C++ con Octave
Octave es un proyecto de software libre con un diseño bastante acertado. Tanto es así que las librerías de extensón de Octave junto con la capacidad de incrustar el intérprete en una aplicación hace que Octave sea una biblioteca de lo más completa. Esto proporciona a Octave una cualidad no tan esperada; puede utilizarse como extensión natural de C++ para el cálculo científico. La cabecera oct.h proporciona todos los tipos necesaros para el cálculo numérico, matrices, vectores fila y columna, operaciones matriciales y vectoriales... ¡Incluso resolución de EDOs! Para abrir boca nada mejor que un pequeño experimento con un programa de lo más simple: /*testoctavestandalone.cc*/
#include <iostream>
#include <oct.h>
int main(void)
{
std::cout << "Octave y c++ \n";
Matrix a=identity_matrix(2,2);
std::cout << a;
}
Para compilarlo utilizaremos la aplicación mkoctfile como viene siendo habitual pero con la opción --link-stand-alone para generar un ejecutable. $ mkoctfile --link-stand-alone testoctavestandalone.cc -o test.exe
Y lo ejecutaremos, en un entorno UNIX... $ ./test.exe
Octave y c++
1 0
0 1
La biblioteca en C++ de Octave es mucho más potente como demuestra el siguiente ejemplo: /*test2.cc*/
#include <iostream>
#include <oct.h>
int main(void)
{
Matrix a=Matrix(2,2);
ColumnVector b=ColumnVector(2);
a(0,0)=2.;
a(1,0)=5.;
a(0,1)=-6.;
a(1,1)=3.;
b(0)=1.;
b(1)=0.;
std::cout << a.solve(b);
return 0;
}
donde resolvemos el siguiente sistema de ecuaciones:
| 漂R> 輯FONT> |
|
?> ?NT> |
x= |
漂R> 輯FONT> |
|
?> ?NT> |
$ mkoctfile --link-stand-alone test2.cc -o test.exe
$ ./test.exe
0.0833333
-0.138889
7.5 Optimización de la evaluación de funciones.(+)
En la sección 5.2.1 hemos planteado la necesidad de optimizar el proceso de evaluación de funciones especialmente complejas. El coste de la evaluación de una función muy compleja puede reducirse significativamente con la interpolación polinómica. Un caso muy llamativo son las funciones con gran cantidad de términos armónicos (senos, cosenos y tangentes) que en la evaluación suelen cancelarse. Aunque la evaluación de funciones no es una de las tareas más optimizadas de Matlab no está de más utilizar los conocimientos adquiridos para aumentar su rendimiento7. Los polinomios de Chebyshev son un recurso bastante habitual porque implican un error de interpolación pequeño. También podríamos pensar en utilizar un desarrollo de Taylor, pero tienen el inconveniente de que sólo aproximan la función en un punto, mientras que lo que nos interesa es un error controlable en todo un intervalo.
Otra característica interesante es que cuando desarrollamos una función en serie, ya sea de Fourier o Polinómica, las operaciones elementales como las derivadas o las integrales son triviales y, lo que es más importante, fáciles de implementar8.
7.5.1 Polinomios de Chebyshev.
Este es un resumen del tratamiento del desarrollo en polinomios de Chebyshev hecho en [3].
Los polinomios de Chebyshev cumplen la siguiente formula:
Tn(x)=cos(n cos-1x)
Esta expresión genera la siguiente serie de polinomios:
T0(x)=1
T1(x)=x
T2(x)=2x2-1
T3(x)=4x3-3x
T4(x)=8x4-8x2+1
⋯
Tn+1(x)=2xTn(x)-Tn-1(x) n≥1
- Estos polinomios tienen la particularidad de contar con n ceros en los puntos:
| x=cos |
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
k=1,2,…,n |
y n+1extremos relativos en:
| x=cos |
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
acotados por el intervalo [-1,1]. Estos polinomios, además, tienen la característica de ser ortogonales según la función peso (1-x2)-1/2. Se demuestra que los coeficientes del desarrollo de la forma:
se calculan mediante:
| cj= |
|
|
f |
漂R> 缂R> 缂R> 輯FONT> |
cos |
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
?> ?> ?> ?NT> |
cos |
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
Evidentemente no necesariamente nos interesa un desarrollo que contenga todos los puntos posibles; lo más normal es calcular un desarrollo truncado de Chebyshev. Con el cambio de variable
conseguimos un desarrollo en el intervalo [a,b], mucho más útil que el [-1,1] original.
La siguiente es una función escrita en C++ que implementa el desarrollo. Es interesante fijarse en cómo escribir la ayuda en formato texinfo /* chebyexpand.cc */
#include <octave/oct.h>
#include <math.h>
static double PI=4.*atan(1.);
DEFUN_DLD(chebyexpand,args, ,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{c}]=} chebyexpand\
(@var{f},@var{a},@var{b},@var{n})\n\
\n\
Chebyshev fit:\n\
\n\
Given a function @var{f} and the interval @var{a},@var{b} computes\n\
the @var{n} coefficients of the Chebyshev polynomial expansion.\n\
\n\
@example\n\
>> chebyexpand(@@(x) 2*x^2-1,-1,1,5)\n\
ans =\n\
-2.6645e-16\n\
1.1948e-17\n\
1.0000e+00\n\
7.9786e-17\n\
-1.0426e-16\n\
@end example\n\
\n\
Note that @var{f} must be a function handle or an anonymous function\n\
\n\
@end deftypefn")
{
int j,k;
octave_value_list tmp;
octave_value_list inval;
octave_function *input_fcn=0;
if (args(0).is_function_handle() || args(0).is_inline_function())
input_fcn = args(0).function_value();
else
{
error("this is not a function handle nor an inline function");
return octave_value(-1);
}
double a=args(1).double_value();
double b=args(2).double_value();
int n=args(3).int_value();
ColumnVector f=ColumnVector(n);
ColumnVector c=ColumnVector(n);
for (k=0;k<n;k++)
{
inval(0)=octave_value(cos(PI*(k+0.5)/n)*\
0.5*(b-a)+0.5*(b+a));
tmp=input_fcn->do_multi_index_op(1,inval);
f(k)=tmp(0).double_value();
}
for (j=0;j<n;j++)
{
double sum=0.;
for (k=0;k<n;k++)
{
sum += f(k)*cos(PI*j*(k+0.5)/n);
}
c(j)=sum*2./n;
}
return octave_value(c);
}
Como ejemplo de cómo debe ser la ayuda de una función de octave iremos al intérprete y teclearemos: -- Loadable Function: [C]= chebyexpand (F,A,B,N)
Chebyshev fit:
Given a function F and the interval A,B computes the N
coefficients of the Chebyshev polynomial expansion.
>> chebyexpand(@(x) 2*x^2-1,-1,1,5)
ans =
-2.6645e-16
1.1948e-17
1.0000e+00
7.9786e-17
-1.0426e-16
Note that F must be a function handle or an anonymous function
Ahora nos falta crear la función que dados los coeficientes nos recupere el valor de la función de un modo eficiente. Para ello utilizamos la siguiente recurrencia también definida en [3]:
dm+2≡ dm+1≡0
dj=2xdj+1-dj+2+cj j=m,m-1,…,2
Pero una función recurrente nos obliga a implementar un bucle, algo que rompe completamente el objetivo de crear una función que sea muy rápida. Sin embargo el esfuerzo de programarla en cualquier otro lenguaje es mínimo9. Por ejemplo, en C++ para octave (crear un archivo .mex sería aún más fácil): El programa tiene un bug #include <octave/oct.h>
DEFUN_DLD(chevyval,args, ,
"Calcula el valor del desarrollo de Chevyshev en x")
{
ColumnVector c (args(0).vector_value());
double a (args(1).double_value());
double b (args(2).double_value());
double x (args(3).double_value());
double xx;
double y;
double ynm1;
double dummy;
int i;
int n =c.length();
xx=(2.*x-a-b)/(b-a);
if ((x-a)*(x-b)>0)
{
error("chevyval: bad input point");
return octave_value_list();
}
else
{
ynm1= c(n-1);
y=2.*xx*ynm1+c(n-2);
for (i=n-3;i<0;i--)
{
dummy=y;
y=2.*xx*y+ynm1+c(i);
ynm1=dummy;
}
y=xx*y-ynm1+0.5*c(0);
return octave_value(y);
}
}
7.6 OCTS (Octave Control Theory Suite)
El control automático es una de las disciplinas más importantes de la ingeniería de sistemas. Se aplica tanto a la aeronautica, el diseño de sistemas mecánicos, el control de redes... Es además una herramienta bastante sencilla con resultados sorprendentes.
A medida que el sistema de control va creciendo se complica su tratamiento analítico pero tratandose de polinomios es sencillo diseñar herramientas numéricas para ello. Cuando afrontamos un problema de control automático solemos empezar con su representación en un diagrama de bloques. Cada uno de los bloques representa una parte de las ecuaciones en forma de funciones de transferencia y encontrar un modo sencillo de crear, manipular y calcular el sistema no es más que un problema de diseño.
Encontramos bastantes programas especializados en el tratamiento de sistemas dinámicos. Podemos dividirlos en dos tipos, las herramientas por línea de comandos y las herramientas gráficas. Las primeras suelen ser tanto o más potentes que las segundas pero su curva de aprendizaje es más larga, incluso puede ser que demasiado larga. En el segundo grupo nos encontramos un programa verdaderamente popular construido sobre Matlab, Simulink.
Simulink fue creado para el tratamiento fácil de sistemas dinámicos pero rápidamente se convirtió en una herramienta de modelado verdaderamente potente. Casi cualquier programa de Matlab puede ser planteado en un diagrama de bloques de Simulink. Diríamos que es la extensión natural de los problemas de control aplicados a Matlab. Simulink es una herramienta fabluosa pero no es parte del lenguaje de programación. Es una aplicación, no una extensión del lenguaje.
Mantener los programas de simulación en el ``espacio del código fuente'' suele ser una práctica bastante beneficisa pero requiere tiempo de desarrollo y experiencia. Este curso pretende ser una ayuda precisamente para adquirir rápidamente esta experiencia. La manipulación de sistemas dinámicos puede ser un buen ejemplo de cómo programar bien en Matlab.
No analizaremos Simulink. Para esto están los libros de Simulink. A cambio trataremos uno de los toolkits más útiles y peor documentados de Octave, La Suite de Teoría del Control de Octave (OCTS).
7.6.1 La representación del sistema
Cada parte de un sistema dinámico puede ser expresada de muchas maneras, en forma de función de transferencia, en el espacio de las variables de estado y proporcionando sus ceros y sus polos.
Una de las particularidades del OCTS es que almacena cada parte del sistema en una estructura de datos. La consecuencia directa de este planteamiento es que para interactuar con lo almacenado estamos obligados a utilizar las funciones diseñadas para ello:
- tf
- Crea un bloque mediante su función de transferencia
- ss
- Crea un bloque mediante su representación en el espacio de sus variabels de estado
- zp
- Crea un bloque mediante la información proporcionada por sus ceros y sus polos
El propio proceso de creación se encarga de calcular las representaciones alternativas y almacenarlas dentro de la estructura. Por ejemplo, para introducir:
| block(s)= |
| K(s+0.5)(s+1) |
|
|
| s2(s+0.1)(s+5)(s+10) | |
con K=1 >> block=zp([-0.5,-1],[-0.1,-5,-10,0,0],1)
El resultado es el siguiente sistema block =
{
inname =
{
[1,1] = u_1
}
k = 1
n = 5
nz = 0
outname =
{
[1,1] = y_1
}
pol =
-0.10000 -5.00000 -10.00000 0.00000 0.00000
stname =
{
[1,1] = x_1
[1,2] = x_2
[1,3] = x_3
[1,4] = x_4
[1,5] = x_5
}
sys =
1 0 1 0
tsam = 0
yd = 0
zer =
-0.50000 -1.00000
}
- sysout
- Muestra las distintas representacoiones del sistema:
Esta función es bastante útil para conocer funciones en tiempo real. Como ya está todo almacenado dentro de la estructura de datos del sistema no tenemos más que solicitar cualquiera de ellas: >> sysout(block,'tf')
Input(s)
1: u_1
Output(s):
1: y_1
transfer function form:
1*s^2 + 1.5*s^1 + 0.5
-----------------------------------------------
1*s^5 + 15.1*s^4 + 51.5*s^3 + 5*s^2 + 0*s^1 + 0
>> sysout(block,'zp')
Input(s)
1: u_1
Output(s):
1: y_1
zero-pole form:
1 (s + 0.5) (s + 1)
------------------------------
s^2 (s + 0.1) (s + 5) (s + 10)
>> sysout(block,'ss')
Input(s)
1: u_1
Output(s):
1: y_1
state-space form:
5 continuous states, 0 discrete states
State(s):
1: x_1
2: x_2
3: x_3
4: x_4
5: x_5
A matrix: 5 x 5
-10.00000 -4.50000 1.00000 0.00000 0.00000
0.00000 -5.00000 1.00000 0.00000 0.00000
0.00000 0.00000 -0.10000 1.00000 0.00000
0.00000 0.00000 0.00000 0.00000 1.00000
0.00000 0.00000 0.00000 0.00000 0.00000
B matrix: 5 x 1
0
0
0
0
1
C matrix: 1 x 5
-9.00000 -4.50000 1.00000 0.00000 0.00000
D matrix: 1 x 1
0
