4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
50
Capítulo 4
Técnicas de reducción de circuitos
4.1 Introducción
La técnica conocida como balance armónico se introdujo anteriormente para resolver las dificultades que presentaban otras técnicas, basadas en transitorios en el dominio del tiempo, al tratar con grandes circuitos no lineales.
La idea principal detrás del balance armónico es que las formas de onda en un circuito no lineal periódicamente excitado, son periódicas, y por tanto, pueden ser representadas en el dominio de la frecuencia como series finitas de Fourier. En efecto, eso transforma las ecuaciones diferenciales que describen el circuito en un conjunto de ecuaciones algebraicas no lineales, que puede resolverse directamente usando técnicas iterativas como el método de Newton-Raphson.
Sin embargo, el principal problema que nos encontramos al resolver estas
ecuaciones no lineales, es el prohibitivo coste computacional requerido para almacenar y factorizar la matriz jacobiana de las ecuaciones no lineales. Desafortunadamente, para circuitos integrados grandes de RF con fuerte no linealidades, esta matriz jacobiana llega a ser enorme y densa, y su factorización puede costarnos mucho incluso para circuitos de un tamaño medio.
Además, como el resto de métodos localmente convergentes, el método de
Newton también tenía problemas de convergencia, como ya vimos. La convergencia puede lograrse sólo encontrando valores iniciales suficientemente cercanos a la solución, lo cuál es difícil, sobre todo en casos de altos niveles de excitación.
Los métodos llamados de continuación, se propusieron para evitar los problemas
de convergencia asociados con el método directo de Newton. En estos métodos, el conjunto de ecuaciones no lineales original es reemplazado por un sistema auxiliar del
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
51
mismo tamaño, cuya solución es trivial, y luego usar algún parámetro de barrido para trazar la trayectoria de la solución de vuelta al sistema inicial. Pero esos métodos con computacionalmente muy costosos, ya que requieren la factorización de la matriz jacobiana para muchos valores del parámetro de barrido. Otras técnicas lineales iterativas, basadas en subespacios de Krylov, también se han usado para tratar de reducir el alto coste de factorizar y almacenar la matriz jacobiana. Sin embargo, esos métodos requieren el uso de condicionadores previos para mejorar su convergencia. Además, aún necesitan operar con el tamaño completo de la matriz jacobiana. Por tanto, en este capítulo se presentan varias formas de aumentar la eficiencia del método Harmonic Newton. Veremos en detalle la estructura del jacobiano y se darán distintas formas de explotar esa estructura. En primer lugar, veremos de forma breve distintas técnicas basadas tanto en métodos exactos como en métodos inexactos de Newton. Posteriormente se mostrará detenidamente cómo es la formación del jacobiano. Y por último, se pasará a ver en mayor detalle algunas de las técnicas introducidas inicialmente. 4.1.1 Métodos de Newton exactos Ya sabemos que la parte más cara del algoritmo, en términos de tiempo de CPU y uso de memoria, es la formulación, inversión y solución de la matriz jacobiana. La inversión del jacobiano es una operación de orden )( 3NO , con N el número de nodos. Hay varias formas de reducir el impacto de este costoso proceso. En primer lugar podríamos citar el método en el que el jacobiano se divide en bloques, donde sólo se usan bloques a lo largo de la diagonal del jacobiano. Fueron Heron, Chang y Steer los que introdujeron este método [9]. Esto permite que sean los bloques individuales los que se inviertan, preferible a invertir la matriz jacobiana completa. Debido a la estructura de la matriz jacobiana, todas las contribuciones lineales están localizadas dentro de estos bloques. También usaron el método de Samanskii, que implica el uso del mismo jacobiano invertido en más de una iteración. Estos son los métodos que se ajustan mejor a sistemas ligeramente no lineales. Fue Rizzoli quien usó técnicas de sparse matrix (matriz formada sobre todo por ceros) para mejorar la simulación. Cuando se almacenan y manipulan en el ordenador matrices de este tipo, a menudo es necesario modificar los algoritmos estándares y aprovecharse de la estructura sparse de la matriz. Por su naturaleza, estas matrices son fácilmente comprimibles, lo que provoca un gran ahorro en el uso de memoria. En las técnicas de Rizzoli, algunos elementos del jacobiano se establecían automáticamente a cero. Se almacenan sólo los elementos distintos de cero, con el índice de su fila y columna. Así se reduce la cantidad de memoria necesaria en calcular el jacobiano, y mejora la velocidad de cálculo. No se realizan operaciones sobre los componentes de la matriz que valen cero, ya que no existen en la memoria. Para tratar simulaciones de circuitos de alta potencia, hay varias técnicas propuestas, entre la que destaca la de reducción de norma, que ajusta la dirección de actualización del método de Newton, de modo que la iteración resultante sea óptima. Rizzoli usó esta técnica en conjunto con el modelado paramétrico de dispositivos no
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
52
lineales, lo cual es útil cuando se emplea la propuesta de las variables de estado. En vez de definir el sistema de incógnitas para ser estrictamente tensiones o corrientes, las incógnitas se eligen que sean variables de estado paramétricas que se mapean en tensiones o corrientes dependiendo de su valor. Otro problema en el contexto de los métodos de Newton y el balance armónico es la simulación de circuitos con un gran número de elementos no lineales. Son los métodos inexactos de Newton los que parecen adaptarse mejor a ese problema. 4.1.2 Métodos de Newton inexactos El principal inconveniente a la hora de usar el método de Newton está en el coste de invertir la matriz del jacobiano. En realidad, el jacobiano no suele invertirse, sino que se hace una descomposición LU para resolver: bJy = (4.1) donde b es el vector de error de la iteración dada. Una vez que se hace la factorización LU, la misma matriz descompuesta puede ser reusada tantas veces como se desee en el método de Samanskii. De todos modos, aunque la inversión de jacobiano no se haga explícitamente, la descomposición del jacobiano sigue siendo un proceso )( 3NO . Para tratar de evitar este cuello de botella computacional, son varios los autores que han propuesto usar técnicas lineales iterativas para resolver (4.1). La ventaja de usarlas es que sólo es necesario multiplicaciones vector-matriz, por lo que el tiempo de resolución aumenta sólo un poco más ligeramente de lo que ocurre en el caso lineal, conforme se incremente el número de incógnitas. Sin embargo, el problema de los métodos lineales es que no convergen fiablemente. Para ayudar en la convergencia, hay que usar algún tipo de precondicionador. Melville fue de los primeros en proponer estas técnicas, y usó un precondicionador (matrix preconditioning) para modificar la ecuación (4.1) a:
bJJxJ 11 ~~ −− = (4.2)
la cual tiene la misma solución. Idealmente, el precondicionador J~
debería ser una buena aproximación a J y también debería ser fácil de invertir. En vez de resolver (4.2), él resolvió:
bzJ =~
(4.3)
Melville también estudió el uso del jacobiano linealizado alrededor del punto de funcionamiento DC del circuito. Esto puede encontrarse fácilmente de la matriz de admitancias del circuito, que debe formarse al principio de la simulación de balance armónico para el cálculo de la respuesta lineal. La matriz jacobiana se organiza de forma que las contribuciones del circuito lineal se localizan en N bloques a lo largo de la diagonal, cuando se usan N frecuencias de análisis. Por tanto, este precondicionador implica inversiones de bloques, lo que es significativamente más cómodo que la
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
53
inversión de la matriz completa, especialmente en circuitos con un gran número de incógnitas. Esta técnica es buena para grandes circuitos que operan en un régimen ligeramente no lineal. Unas de las familias de técnicas más destacadas son las de subespacios de Krylov. Estos métodos se caracterizan por los subespacios en los que se encuentran las iteraciones. El subespacio de Krylov k-ésimo para una matriz dada J y un vector b es bJbJJbb k 12 ,...,,, − . Si empezamos con una estimación )0(x proporcional a b, la
iteración produce aproximaciones )( jx que son combinaciones de los vectores del subespacio. Hay distintas técnicas basadas en estos subespacios como GMRES (residuo mínimo generalizado), MINRES (residuo mínimo), CG (gradiente conjugado),... Feldman y Melville desarrollaron algoritmos basados en subespacios de Krylov, y para tratar de manejar circuitos no lineales más generales, hicieron mejoras en las matrices de precondicionamiento. En vez de usar sólo las contribuciones lineales del jacobiano, se usaban también las contribuciones no lineales que aparecían en los bloques de la diagonal del jacobiano. Aunque lo cierto es que para fuertes no linealidades, también debe incluirse en el precondicionador alguna información de bloques de fuera de la diagonal. Rizzoli también exploró el uso de métodos inexactos de Newton para grandes circuitos de microondas. Él compara los métodos exactos e inexactos escribiendo: iii nXX +=+1 (4.4)
donde in es la actualización del método exacto. Para un gran número de incógnitas,
Rizzoli propuso que puede ser más eficiente usar un método inexacto, con una actualización is , de modo que:
)()()( iiiii XEfsXJXE ≤+ (4.5)
y iii sXX +=+1 (4.6)
donde 10 <≤ if , cuando 0=if , la actualización se reduce a la actualización exacta. Y
en otro caso, if sirve como una medida de cuánto difiere la actualización inexacta de la
exacta. La elección de if se va actualizando en cada iteración. Se empieza con 5.00 =f
y se calculan los siguientes mediante:
)(
)()()(
1
111
−
−−− −−=
i
iiii
iXE
sXJXEXEf (4.7)
Una vez se elige if , se resuelve aproximadamente la ecuación de Newton hasta que se
satisface (4.4). Es la técnica GMRES (residuo mínimo generalizado) la que se usa para calcular la solución aproximada. Y para mejorar las propiedades de convergencia
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
54
también hay que seleccionar un precondicionador adecuado. Rizzoli reemplaza la ecuación exacta de Newton: )()( iii XEnXJ −= (4.8)
con )()( 1
iiiii XEnPPXJ −=− (4.9)
Es importante que P sea una aproximación cercana a J y que sea fácil de invertir. Una vez que se define un precondicionador, la estimación inicial para )0(
is se define por:
)(1)0(
iii XEPs −−= (4.10)
y se define un conjunto de vectores reales de longitud N: [ ] )()( 1)1(
iiiNi XEPXJIK −−= (4.11)
11)( )( −−= q
iii
q
i KPXJK (4.12)
donde NI es la matriz identidad. El espacio vectorial extendido por los vectores )(q
iK ,
Qq ≤≤1 se llama subespacio de Krylov de dimensión Q. La aproximación de orden Q
de is viene dada por:
∑=
−−+=Q
q
q
iqii
Q
i KPss1
11)0()( α (4.13)
Si la ecuación (4.4) se satisface, )(Q
is se toma como la actualización aproximada
de Newton (lo cual está garantizado para Q suficientemente grande, ya que
i
Q
iQ ns =∞→)(lim ). Rizzoli propone usar una dimensión del espacio de Krylov más
pequeña en las primeras actualizaciones, e irla incrementando gradualmente conforme se alcanza la solución final, para lo que 50≤Q es suficiente para la mayoría de problemas de balance armónico. 4.1.3 Formación de la matriz jacobiana Como se discutió en el Capítulo 3, lo primero que hay que hacer con el circuito no lineal es dividirlo en un subcircuito lineal y otro no lineal. Esto resulta en un subcircuito no lineal de N puertos, cada uno de ellos conectado a un elemento no lineal. La solución DC suele encontrarse primero, y usarse luego como la estimación inicial. Ya que no se requieren cálculos de FFT y que el número de variables se reduce ampliamente, la solución DC es generalmente fácil de encontrar a través del método de Newton.
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
55
Al final, de lo que se trataba era de evaluar la expresión (3.34), que representa una iteración del método de Newton. Para ello, es necesario calcular )(VF y )(VJ F . Expresemos la ecuación de balance armónico de este modo:
0)~
(~
)~
(~
)( 11 =+Ω++= −− VFiFVFqFjYVIVF GS (4.14)
donde F~
es el símbolo de la transformada de Fourier. Este sistema no lineal de ecuaciones se resuelve tal que ε<)(VF , para algún ε suficientemente pequeño
definido por el usuario. En F(V), la respuesta lineal se calcula directamente mediante multiplicaciones de matrices y vectores. Y para evaluar los términos no lineales en )(VF , los espectros (dominio de la frecuencia) de las tensiones de los nodos son transformados al dominio del tiempo (aplicando FFT inversa), aplicados luego a los dispositivos no lineales, y las resultantes formas de onda (dominio del tiempo) de las corrientes se transforman al dominio de la frecuencia (mediante FFT). En cambio, el cálculo de la matriz jacobiana )(VJ F es más complicado. El jacobiano consiste en derivadas de la función de error F(V) con respecto a las tensiones incógnitas. Las derivadas de la parte lineal de F(V) se encuentra fácilmente, y consiste en la matriz de admitancias lineal. Pero las contribuciones de la parte no lineal no son tan fáciles de calcular. En parte porque )(vf y v están obligadas a ser funciones reales (formas de onda en el dominio del tiempo tienen que serlo), y también porque deseamos no usar frecuencias negativas. Aplicar estas restricciones resulta en que )(VJ F no sea representable en el campo complejo. Para solucionar este problema, cada número complejo se escribe como un vector equivalente en 2ℜ . De modo que si tenemos un número complejo X, podemos definir números reales RX y IX ( )Re(XX R = ,
)Im(XX I = ) y un vector 2ℜ∈X tal que [ ]TIR XXX ,= . Una notación similar a esta es la que usaremos para los vectores, funciones y matrices. Según lo que acabamos de decir, vamos a evitar usar las frecuencias negativas, y gracias a eso, el jacobiano pasará de ser una matriz compleja a una matriz real, aunque sea de dimensión mayor. Esta conversión reduce la memoria requerida y el número de operaciones necesarias. Ahora el jacobiano se define así:
V
VI
V
VQY
V
VFVJ G
NxNF∂
∂+
∂
∂Ω+=
∂
∂=
)()()()( (4.15)
∂
∂
∂
∂
∂
∂
∂
∂
=∂
∂=
)(
),(
)(
),()(
),(
)(
),(
)(
),(),,(,
lV
kVF
lV
kVF
lV
kVF
lV
kVF
lV
kVFlkVJ
I
n
I
R
n
I
m
I
n
R
m
R
n
R
m
n
m
mnF (4.16)
−=
),(),(
),(),(),(
lkYlkY
lkYlkYlkY
R
mn
I
mn
I
mn
R
mn
mn (4.17)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
56
para Nnm ,....,2,1, ∈ .
−=Ω
0
0
0
0
kw
kw (4.18)
Así, las componentes real e imaginaria de las variables complejas son almacenadas de forma separada como números reales, requiriendo el cálculo de cuatro derivadas para describir la relación entre corriente y tensión o carga y tensión a cualquier frecuencia dada.
Por tanto, para el cálculo de V
VIG
∂
∂ )( y
V
VQ
∂
∂ )( no nos valen las expresiones que
vimos en (3.49). Hay que realizar cuatro derivadas para calcular cada una de ellas.
Vamos a ver cómo sería )(
),(,
lV
kVIR
n
R
mG
∂
∂. El resto de términos en
V
VIG
∂
∂ )( y los de
V
VQ
∂
∂ )(
se harían de forma similar. Partiendo de la definición de transformada de Fourier se llega a:
∫=T
PmgmG dttktiT
kVI0 ,, )cos()(
1),( ω (4.19)
Al aplicar la regla de la cadena:
∫ ∂
∂
∂
∂=
∂
∂ T
P
n
n
n
mg
R
n
R
mGdttk
lV
tv
tv
ti
TlV
kVI
0
,, )cos()(
)(
)(
)(1
)(
),(ω (4.20)
Ahora expresamos )(tvn en términos de parte real e imaginaria:
∑∞
=
−+=1
)()()cos()(2)0()(k
P
I
nP
R
n
R
nn tksenkVtkkVVtv ωω (4.21)
Y teniendo en cuenta que para 0≠l :
−=
∂
∂∂
∂
=∂
∂
)(2
)cos(2
)(
)()(
)(
)(
)(
tlwsen
tlw
lV
tv
lV
tv
lV
tv
P
P
I
n
n
R
n
n
n
n (4.22)
y para 0=l :
=
∂
∂
0
1
)0(
)(
n
n
V
tv (4.23)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
57
Para 0≠l queda:
=∂
∂=
∂
∂∫T
PP
n
mg
R
n
R
mGdttltk
tv
ti
TlV
kVI
0
,, )cos()cos()(
)(2
)(
),(ωω
[ ]∫ −++∂
∂=
T
PP
n
mgdttlktlk
tv
ti
T 0
, ))cos(())cos(()(
)(1ωω (4.24)
y para 0=l :
∫ ∂
∂=
∂
∂ T
P
n
mg
R
n
R
mGdttk
tv
ti
TV
kVI
0
,, )cos()(
)(1
)0(
),(ω (4.25)
Y teniendo en cuenta la definición de )(kGmn que vimos en (3.47) y el resto de
derivadas, tendré, para 0≠l :
+−−++−
−−+++−=
∂
∂
)()()()(
)()()()(
)(
),(,
lkGlkGlkGlkG
lkGlkGlkGlkG
lV
kVIR
mn
R
mn
I
mn
I
mn
I
mn
I
mn
R
mn
R
mn
n
mG (4.26)
y para 0=l :
=
∂
∂
0)(
0)(
)0(
),(,
kG
kG
V
kVII
mn
R
mn
n
mG (4.27)
La derivada V
VIG
∂
∂ )( será una matriz formada por bloques como esos. Y para
V
VQ
∂
∂ )( es similar pero con C en vez de G.
Hay que destacar que es preciso forzar que 0)Im( =dcV , ya que los fasores
deben ser reales en continua. Como las cantidades DC son estrictamente reales, no requieren tanto espacio de almacenamiento en el jacobiano.
Podría ocurrir que en estas matrices de conversión hubiese filas y columnas consistentes sólo en ceros (segunda fila y segunda columna). Esto conlleva a que haya una singularidad en la estructura del jacobiano, y por tanto, al ser singular, no tiene inversa, y no puede aplicarse el método de Newton. Si esto ocurriese, podríamos deshacernos de esta singularidad eliminando esa fila o columna, o bien haciendo que en la matriz original (antes de la conversión), la diagonal principal no tenga ceros [1].
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
58
4.2 Técnicas de aproximaciones del jacobiano y
precondicionadores
Con lo visto en el apartado anterior, y como se ve en (4.16), el jacobiano es organizado como una matriz de bloques. Cada bloque es una matriz de conversión que a su vez es una matriz de bloques, consistente en bloques 2x2 que resultan de que los coeficientes de Fourier pertenezcan a 2ℜ en vez de al plano complejo. Las matrices de conversión son completas si están asociadas con un nodo que tiene conectado un dispositivo no lineal. En caso contrario son diagonales. En un circuito MMIC, hay dispositivos no lineales conectados a la mayoría de los nodos, por lo que las matrices de conversión suelen estar completas. Claramente, el jacobiano armónico puede ser muy grande y mucho más denso que la mayoría de matrices de circuitos. Aplicar técnicas tradicionales de sparse matrix no es suficiente para resolver eficientemente la ecuación de actualización del método de Newton. Es necesario también reducir la densidad de la matriz. El jacobiano se usa sólo para generar nuevas iteraciones; no se usa para confirmar la convergencia, por lo que los errores en el jacobiano afectan sólo a la tasa y a la región de convergencia, no a la precisión de la solución final. Las aproximaciones en el jacobiano reducen esa tasa de convergencia, pero lo que se gana en eficiencia es más importante. Quizás, la aproximación más sencilla que puede hacerse es simplemente reusar el jacobiano de una iteración anterior. Si el circuito está comportándose casi linealmente, el jacobiano no variará mucho de una iteración a otra. Esta idea, que se atribuye a Samanskii, puede reducir ampliamente el tiempo requerido para una iteración, porque elimina completamente la construcción y descomposición LU del jacobiano. Por tanto, sólo se necesitan los pasos de sustitución hacia delante y hacia atrás. La segunda aproximación del jacobiano resulta de explotar las características naturales de las matrices de conversión de los dispositivos no lineales. Sabemos que los algoritmos usuales de álgebra lineal para matrices generales, requieren una cantidad de operaciones de orden cúbico en el tamaño de la matriz en consideración. Sin embargo, en el caso en que dicha matriz sea estructurada en algún sentido, la complejidad de su resolución se reduce de forma drástica. Así, las matrices de conversión pueden dividirse en la suma de una matriz Toeplitz (con la forma dada por jiij ra −= ) y una matriz
Hankel ( jiij sa += ). Por ejemplo:
SRV
VI
n
m +=∂
∂ )( (4.28)
donde
=−
−−
−−−
0123
1012
2101
3210
rrrr
rrrr
rrrr
rrrr
R
=
6543
5432
4321
3210
ssss
ssss
ssss
ssss
S
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
59
−=
)()(
)()(
kGkG
kGkGr
R
mn
I
mn
I
mn
R
mn
k
−=
)()(
)()(
kGkG
kGkGs
R
mn
I
mn
I
mn
R
mn
k (4.29)
donde mnG es la suma de las derivadas de los espectros de resistores no lineales entre el
nodo m y n, o entre m y tierra si nm = . Este espectro tiene la característica de que cuanto más lineales se están comportando los dispositivos que lo generan, más domina la componente DC al resto de armónicos, y más rápido decae la magnitud a altos armónicos. Como resultado, los elementos en la matriz de conversión lejos de la diagonal, serán pequeños comparados a los de la diagonal. Para reducir la densidad del jacobiano, estos pequeños términos lejanos a la diagonal se ignorarán [1].
Varios autores, mencionados en la introducción, han usado precondicionadores para aproximar el jacobiano en técnicas lineales iterativas. La elección del precondicionador o de la aproximación del jacobiano equivale a la selección de qué elementos del jacobiano eliminar antes de la descomposición de la matriz. Hay varias opciones implementadas, pero todas ellas afectan sólo a las contribuciones no lineales del jacobiano, ya que las contribuciones lineales son necesarias para la convergencia, y no se eliminan en ninguna de las técnicas de precondicionamiento.
Una técnica sería usar sólo las contribuciones lineales del jacobiano. Ya que las
contribuciones lineales son constantes con respecto a los valores de las incógnitas, el jacobiano se calcularía y descompondría sólo una vez, al principio de la simulación. La contribución no lineal en el jacobiano no tendría que calcularse, y cómo se calculan mediante la FFT, eliminar la necesidad de este cálculo mejora el tiempo de simulación. La matriz descompuesta puede reusarse en cada iteración del método de Newton siguiendo la idea de Samanskii. Además, debido a que toda la información del jacobiano se almacena en bloques a lo largo de la diagonal, los bloques pueden descomponerse individualmente, ahorrando tiempo. Esta técnica se implementa fácilmente, y puede funcionar en sistemas débilmente no lineales, pero no es robusta. Circuitos con no linealidades más fuertes, requieren alguna información no lineal en el jacobiano para lograr la convergencia del método de Newton.
Otra técnica de aproximar el jacobiano sería usar sólo las contribuciones no
lineales que ocurren en la diagonal de la matriz. En este método no se consideran interacciones entre cantidades del circuito a distintas frecuencias ni a distintos nodos (por ejemplo, la relación entre dI y gsv en un FET no se considera, ya que no
pertenecen al mismo nodo). Añadir sólo contribuciones no lineales de la diagonal al jacobiano lineal, no añade complejidad al proceso de descomponer la matriz, ya que todos los elementos de la diagonal contendrán ya contribuciones del circuito lineal. Lo negativo es que debido a los elementos no lineales, el jacobiano tendría que volver a ser recalculado. Lo positivo es que como la información del jacobiano sigue estando en bloques en la diagonal, esto permite una rápida descomposición LU. Así, mientras que este método es más robusto que no usar ninguna contribución no lineal, aún no es suficiente para circuitos fuertemente no lineales, ya que no se está considerando información entre distintas frecuencias.
El método llamado Block Jacobian expande el rango de contribuciones no
lineales consideradas en el jacobiano. El tamaño del bloque puede variarse de modo que
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
60
más relaciones entre frecuencias sean consideradas, permitiendo una representación más precisa del jacobiano, pero manteniendo la estructura de bloque, que reduce el tiempo de descomposición LU. En la estructura del jacobiano, las contribuciones de elementos lineales ocurren sólo entre variables a la misma frecuencia. Sólo los elementos no lineales pueden producir corrientes a frecuencias distintas de las de las tensiones de excitación. Si el jacobiano se estructura de modo que las relaciones entre los parámetros del circuito sean agrupadas por frecuencias, todas las contribuciones lineales del jacobiano estarán en bloques a lo largo de la matriz diagonal. Los únicos bloques fuera de la diagonal serán los de contribuciones no lineales. Por ejemplo, un circuito con N nodos y K frecuencias de análisis tendrá una matriz con la siguiente estructura:
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
=
)(
),(
)0(
),(
)0(
),(
)0(
)0,(
)0(
)0,(
)0(
)0,(
)0(
)0,(
)0(
)0,(
)(
)0,(
)0(
)0,(
)0(
)0,(
)0(
)0,(
)(
11
11
1
1
1
1
11
1
1
1
1
KV
KVF
V
KVF
V
KVF
V
VF
V
VF
V
VF
V
VF
V
VF
KV
VF
V
VF
V
VF
V
VF
VJ
I
N
I
N
I
I
N
R
I
N
I
N
I
N
I
I
N
R
I
N
I
I
R
I
I
N
R
I
N
R
I
R
R
R
F
LLL
MMMLMM
MML
MMMMMM
LLLL
LL
(4.30)
Como se explicó antes, la matriz de admitancias lineal contribuye sólo a los )(
),(
hV
iVF
n
m
∂
∂
donde hi = . Esos elementos ocurren sólo en bloques a lo largo de la diagonal de la matriz. Podemos rescribir la matriz en términos de bloques de distinta frecuencia:
=
KNNNN
K
K
K
F
BBBB
BBBB
BBBB
BBBB
VJ
,2,1,0,
,22,21,20,2
,12,11,10,1
,02,01,00,0
)(
L
MLMMM
L
L
L
(4.31)
donde cada bloque ijB contiene todas las derivadas de las funciones a la frecuencia i-
ésima con respecto a las variables de estado a la frecuencia j-ésima. Estas técnicas de jacobiano en bloques pueden dividirse en dos subcategorías, de bloques en la diagonal y de bloques fuera de la diagonal. La primera es cuando sólo se
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
61
usan los bloques de la diagonal. Los elementos usados de cada columna y cada fila representan sólo un bloque. Vemos dos ejemplos en las Figura 4.1.
Figura 4.1: a) 16 bloques en la diagonal; b) 4 bloques en la diagonal. Cada bloque pequeño en b) representa todas las derivadas de una función a una frecuencia respecto a todas las variables de estado a otra frecuencia dada. En este método de bloques en la diagonal, el número de bloques usados puede variar desde dos hasta el número de frecuencias de análisis, incrementando en potencias de dos. Cada bloque contiene información relacionada con el mismo número de componentes de frecuencia. Por ejemplo, si hay 8 frecuencias de análisis, podemos tener 2, 4 u 8 bloques disponibles. Si se eligen 8, cada bloque tendrá información referente a una sólo frecuencia, mientras que si se eligen 4 bloques, cada uno tendrá información sobre 2 frecuencias. Las técnicas con bloques fuera de la diagonal se basan en la proximidad de un bloque dado a los bloques de la diagonal cuando se está usando el mayor nivel de bloques en la diagonal. En ese caso, hay K bloques diagonales para un circuito con K frecuencias de análisis. Si los índices de la fila y columna del bloque son i y j respectivamente, una técnica con nivel q de bloques fuera de la diagonal incluiría todos los bloques para los que qji ≤− . En la Figura 4.2 se ven algunos ejemplos.
Figura 4.2: a) Nivel 1 de bloques fuera de la diagonal; b) Nivel 3.
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
62
Para circuitos con no linealidades moderadas, el método del jacobiano con bloques en la diagonal supera, generalmente, al resto de métodos. Cuando el número de bloques es igual al número de frecuencias de análisis, las contribuciones extra no lineales suponen una pequeña diferencia en el tiempo total de ejecución y en requerimientos de almacenamiento de memoria comparado con la matriz lineal. Cuanto mayores sean las no linealidades, mayores deberán ser los bloques para que se obtenga convergencia, y el incremento en el tiempo necesario para la ejecución es proporcional al incremento del tamaño de bloques. Usar las técnicas de bloques fuera de la diagonal en estas situaciones, puede ayudar a la obtención de convergencia con una mejora en el tiempo necesario con respecto a cuando se usa la matriz jacobiana entera, aunque la pérdida de la estructura con bloques en la diagonal limita esa mejora [8].
4.3 Técnica basada en subespacios de Krylov 4.3.1 Introducción
Varios algoritmos son los que se han propuesto en la literatura técnica reciente para hacer frente a algunos aspectos de los problemas de simulación. Métodos híbridos tiempo-frecuencia se han introducido para manejar transitorios y señales digitalmente moduladas. Con estas técnicas, el análisis de régimen permanente en microondas se desacopla del análisis de envolvente (banda base), de modo que la simulación se reduce a la solución separada de una secuencia de problemas de multitono, cada uno asociado con uno de los instantes de muestreo de la envolvente.
Los problemas de análisis multitono con un gran número de dispositivos no
lineales, han sido enfrentados eficientemente con algoritmos de balance armónico tales como GMRES (residuo mínimo generalizado). Esta técnica suprime la necesidad de almacenar y factorizar la matriz jacobiana, y reduce la mayor parte del proceso de resolver el sistema a una secuencia de multiplicaciones matriz-vector, que puede ejecutarse eficientemente con la FFT. Métodos de Newton inexactos y de balance armónico a trozos también se han combinado para minimizar eficientemente el número de iteraciones y el de incógnitas.
Para aplicaciones CAD de propósito general, es probable que la explotación óptima de los recursos disponibles del ordenador se logre con la integración de análisis orientado a envolvente con métodos de subespacios de Krylov. Es eso lo que se va a pretender, un algoritmo para la aplicación de un método de Newton inexacto basado en GMRES a un análisis orientado a envolvente de grandes subsistemas de microondas no lineales. Esta técnica fue propuesta por Vittorio Rizzoli et al. [10] 4.3.2 La técnica MIHB
Vamos a ver un algoritmo que realiza análisis de balance armónico orientado a modulación, que usa métodos inexactos de Newton para resolver el sistema no lineal. Esta aproximación se llamará, por brevedad, técnica MIHB (modulation-oriented
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
63
inexact Newton harmonic-balance), y está especialmente aconsejada para la simulación de grandes circuitos no lineales forzados por portadoras de RF/microondas lentamente moduladas. Formulación del problema
Vamos a asumir que las tensiones de las fuentes son señales moduladas en amplitud y de fase cuasiperiódica:
∑ Ω=k
tj
kketFtf )()( (4.32)
donde kΩ es una intermodulación genérica, producto de un conjunto de frecuencias
fundamentales hw ( Hh ≤≤1 ):
∑=
=ΩH
h
hhk wk1
(4.33)
Las no linealidades del circuito transferirán la modulación a todas las formas de
onda de señal soportadas por el circuito, y en particular a las variables de estado. El vector de estado x(t) tendrá la siguiente expresión:
∑ Ω=k
tj
kketXtx )()( (4.34)
La cantidad compleja )(tX k se llamará modulación compleja (o envolvente), o
fasor dependiente del tiempo del producto de intermodulación k-ésimo. También será referido como el k-ésimo armónico de x(t) dependiente del tiempo. El conjunto de todos los armónicos )(tX k se llamará espectro de x(t).
Asumiremos que )(tFk , )(tX k son funciones arbitrarias limitadas en banda, para
los que la representación integral de Fourier siempre existe:
∫Ω
Ω−=
B
B
dwewFtF jwt
kk )()( (4.35)
∫Ω
Ω−=
B
B
dwewXtX jwt
kk )()( (4.36)
donde )(wFk , )(wX k no son dependientes del tiempo. El conjunto de todas las
cantidades complejas )(wX k representa el espectro físico de x(t). BΩ representa el
mayor de los anchos de banda de modulación de todas las señales de interés. El análisis de la técnica MIHB a desarrollar se sostiene sobre la base de que las funciones de modulación son banda base. Es decir, el límite de la banda por arriba, BΩ , debe ser pequeño con respecto a los productos de intermodulación de las RF fundamentales a considerar en el análisis:
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
64
kB Ω<<Ω (4.37)
En relación a esta última expresión, hay que remarcar un par de cosas. La
primera es que en algunos casos, el algoritmo continuará funcionando, en el sentido numérico, incluso cuando esa expresión no sea estrictamente cierta para todo k, pero el resultado podría llegar a ser cuestionable desde el punto de vista físico. La segunda, es que esa ecuación no se cumple, lógicamente, para k=0 ( 00 =Ω ). Esto implica que la
componente DC requiere un tratamiento especial, como veremos más adelante. MIHB usa los armónicos dependientes del tiempo )(tX k como las incógnitas
del problema. Siendo más precisos, )(tX k se muestrean a un número de instantes de
tiempo uniformemente espaciados nt ( Nn ≤≤1 ), y las cantidades complejas )( nk tX
son las elegidas como incógnitas. Los instantes de tiempo de muestreo deben elegirse de modo que satisfagan el teorema de muestreo:
B
tΩ
<∆π
(4.38)
Asumiremos que la información de entrada consiste exclusivamente en los
valores (complejos) )( nk tF .
Ecuaciones de las subredes lineal y no lineal
De acuerdo a la técnica de balance armónico a trozos, el circuito no lineal se divide en una subred lineal y otra no lineal interconectadas a través de Dn puertos, como se muestra en la Figura 4.3. Por simplicidad, se asume que la subred no lineal puede ser descrita por el siguiente conjunto de ecuaciones paramétricas:
= )(,
)(,
)(),()(
2
2
txdt
txd
dt
tdxtxutv d , (4.39)
= )(,
)(,
)(),()(
2
2
txdt
txd
dt
tdxtxuti d , (4.40)
donde v(t) e i(t) son vectores de tensiones y corrientes de los puertos, y )(txd es un
vector de variables de estado de retrasos de tiempo ( )()( mmdm txtx τ−= , con mτ una
constante). Todos los vectores de esas expresiones tienen tamaño Dn .
A una frecuencia angular genérica Ω , la subred lineal puede describirse por una ecuación en el dominio de la frecuencia de la siguiente forma:
)()()()()( ΩΩ+ΩΩ=Ω FYVYI TL (4.41)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
65
donde )(ΩY es la matriz de admitancias en los puertos del dispositivo, cuando todos los
puertos externos están cortocircuitados, )(ΩTY es la matriz de transadmitancias de la subred lineal ampliada (incluyendo resistencias de fuente y carga) desde los puertos externos hasta los del dispositivo, y )(ΩF es el vector de fasores complejos de las fuentes de tensión sinusoidal de frecuencia angular Ω , conectadas a puertos externos.
Si no hay fuentes conectadas a un puerto externo dado, la entrada
correspondiente de )(ΩF es cero, y el puerto está cortocircuitado. Las corrientes entrando en los puertos del circuito de la subred lineal tendrán el subíndice L, y se llamarán convenientemente corrientes lineales. )(ΩLI es así un vector Dn de fasores de corriente lineal a la frecuencia Ω .
Figura 4.3: Representación esquemática de la subred lineal ampliada
El análisis MIHB debe resolver simultáneamente las ecuaciones (3.39), (4.40) y (4.41) de las subredes no lineal y lineal, bajo un conjunto de excitaciones de la forma indicada por (4.32).
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
66
Cálculo de la respuesta de la subred no lineal
Para encontrar la respuesta de x(t) de la subred no lineal, MIHB hace uso de una aproximación cuasiestacionaria. La expresión (4.37) implica que )(tFk , )(tX k son muy
lentas en comparación a las señales de microondas RF. Las constantes de tiempo del circuito RF son tales que, para cada n, puede alcanzarse un régimen permanente en un tiempo mucho más corto que t∆ . Así, el régimen eléctrico del circuito bajo una modulación RF puede describirse como una secuencia de regímenes permanentes de microondas RF, cada uno asociado a un conjunto de valores que las envolventes y sus derivadas toman en uno de los instantes de tiempo de muestreo. Para desarrollar este concepto cuantitativamente, hay que reorganizar las expresiones (4.39) y (4.40):
∑ Ω
+Ω=
k
tjkkk
kedt
tdXtXj
dt
tdx )()(
)( (4.42)
∑ Ω
+Ω+Ω−=
k
tjkk
kkkke
dt
tXd
dt
tdXjtX
dt
txd2
22
2
2 )()(2)(
)( (4.43)
Así, el vector de variables de estado de retrasos de tiempo puede escribirse:
tjjwt
k
BjwT
k
Tj
dk
B
k edwewXeetxΩ
Ω
Ω−
−Ω− ⋅
⋅= ∫∑ )()( (4.44)
donde T es la matriz diagonal de los retrasos de tiempo mτ . Los retrasos de fase
asociados con las constantes de tiempo mτ serán normalmente muy pequeñas en banda
base ( 1<<Ω mBτ ), por lo que para una frecuencia genérica w banda base ( Bw Ω≤ ),
tenemos:
jwTe jwT −≈− 1 (4.45) Las ecuaciones anteriores quedan:
tjkk
k
Tj
dkk e
dt
tdXTtXetx
ΩΩ− ⋅
−≈∑
)()()( (4.46)
Para calcular (4.42), (4.43) y (4.46), hay que conocer las derivadas dt
tdX k )(y
2
2 )(
dt
tXd k . Estas derivadas pueden ser evaluadas en un instante de tiempo de muestro nt ,
de la siguiente forma:
∑=
−
=
≈M
m
mnkm
tt
k tXadt
tdX
n0
)()(
(4.47)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
67
∑=
−
=
≈M
m
mnkm
tt
k tXbdt
tXd
n0
2
2
)()(
(4.48)
donde los coeficientes ma y mb pueden encontrase explícitamente en muchos libros
matemáticos.
En principio, la precisión de las derivadas y del análisis, debería incrementar para un número dado de puntos de muestreo. En la práctica, no se observan cambios apreciables en los resultados para M>2. Como el tiempo consumido por la CPU es una función que incrementará con M, M=3 será el valor escogido por defecto para el análisis MIHB.
La respuesta de la subred no lineal puede calcularse ahora con una DFT multidimensional. Introduciremos las variables de tiempo normalizadas:
wtz = ,
twz hh = ( Hh ≤≤1 ) (4.49)
Obtenemos:
∑∑
= =
k
zkj
k
H
h
hh
etXtx 1)()( (4.50)
Por tanto, queda:
∫Ω
Ω−=
B
B
dwewXtX jz
kk )()( (4.51)
Debido a (4.37), z cambia muy lentamente con respecto a cada hz , y tiene que
asumirse que permanece virtualmente constante en cualquier intervalo de longitud
hw/2π ( Hh ≤≤1 ). Así, alrededor de un instante de muestreo genérico nt , en las
ecuaciones de la subred no lineal ponemos las expresiones que dependen del tiempo a través sólo de las hz :
∑∑
=→ =
k
zkj
nkn
H
h
hh
etXtxtx 1)()()( (4.52)
∑ ∑∑
+Ω=→ =
=−
k
zkjM
m
mnkmnkkn
H
h
hh
etXatXjtydt
tdx1
01 )()()(
)( (4.53)
∑ ∑∑∑
+Ω+Ω−=→ =
=−
=−
k
zkjM
m
mnkm
M
m
mnkmknkkn
H
h
hh
etXbtXajtXtydt
txd1
00
222
2
)()(2)()()(
(4.54)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
68
∑ ∑∑
⋅
−=→ =
=
−
Ω−
k
zkjM
m
mnkmnk
Tj
nd
H
h
hh
k etXaTtXetdtx 1
0
)()()()( (4.55)
Y ahora, expresemos las tensiones y corrientes de los puertos de la subred no
lineal, v(t) e i(t), como series de Fourier, introduciendo los conceptos de espectros dependientes del tiempo, )(tU k y )(tWk :
∑ Ω=k
tj
kketUtv )()( (4.56)
∑ Ω=k
tj
kketWti )()( (4.57)
Y alrededor de un instante de muestreo genérico nt , tenemos:
∑∑
≈ =
k
zkj
nk
H
h
hh
etUtv 1)()( (4.58)
∑∑
≈ =
k
zkj
nk
H
h
hh
etWti 1)()( (4.59)
Como vemos en esas expresiones, la dependencia con el tiempo es sólo a través
de hz , así que )( nk tU y )( nk tW pueden ser calculadas con DFTs de dimensión H
(basado en algoritmo FFT) con respecto a únicamente hz . Para hacer esto, primero las
variables hz son muestreadas uniformemente entre [0 2π ], y los valores muestreados
de )(txn , )(1 ty n , )(2 ty n , y )(td n se calculan mediante X. La FFT multidimensional con
respecto a las hz da directamente )( nk tU y )( nk tW .
Cálculo de la respuesta en RF de la subred lineal
El vector de corrientes lineales entrando en los puertos puede escribirse así:
∑ Ω=k
tj
LkLketIti )()( (4.60)
De la ecuación (4.41) de la red lineal, obtenemos:
[ ]∫Ω
Ω−+Ω++Ω=
B
B
wFwYwUwYtI kkTkkLk )()()()()( dwe jwt (4.61)
Esta expresión no parece útil, porque el espectro físico )(wU k no es conocido.
Sin embargo, bajo las mismas presunciones del apartado anterior, si w es una frecuencia en banda base que satisface Bw Ω≤ , podemos hacer las siguientes aproximaciones
para todo k distinto de cero:
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
69
2
2
2 )()()()( w
d
Ydw
d
dYYwY
kk
kk
Ω=ΩΩ=Ω Ω
Ω+
Ω
Ω+Ω≈+Ω (4.62)
2
2
2 )()()()( w
d
Ydw
d
dYYwY
kk
TTkTkT
Ω=ΩΩ=Ω Ω
Ω+
Ω
Ω+Ω≈+Ω (4.63)
Sustituyendo:
−Ω+Ω
Ω−
Ω
Ω−Ω≈
Ω=ΩΩ=Ω
)()()()()()(
)()()(2
2
2
2
tFYdt
tUd
d
Yd
dt
tdU
d
dYjtUYtI kkT
kk
kkLk
kk
2
2
2
2 )()()()(
dt
tFd
d
Yd
dt
tdF
d
dY kTkT
kk Ω=ΩΩ=Ω Ω
Ω−
Ω
Ω− (4.64)
Y si lo calculamos para el instante de muestreo nt , y haciendo uso de las
aproximaciones de las derivadas que vimos anteriormente, queda:
+⋅Ω
Ω−⋅
Ω
Ω−Ω≈ ∑∑
=−
Ω=Ω=−
Ω=Ω
M
m
mnkm
M
m
mnkmnkknLk tUbd
YdtUa
d
dYjtUYtI
kk0
2
2
0
)()(
)()(
)()()(
∑∑=
−
Ω=Ω=−
Ω=Ω
⋅Ω
Ω−⋅
Ω
Ω−Ω+
M
m
mnkmT
M
m
mnkmT
nkkT tFbd
YdtFa
d
dYtFY
kk0
2
2
0
)()(
)()(
)()(
(4.65) Cálculo de la respuesta en DC de la subred lineal
La expresión (4.65) funciona con suficiente precisión para todo 0≠k . Para 0=k la situación puede ser bastante diferente, porque las constantes de tiempo del
circuito DC son normalmente grandes, en comparación a las del circuito RF. Esto puede hacer que alguna de las aproximaciones anteriores sea errónea. Para solucionar esto, la componente DC de la respuesta de la subred lineal se calcula mediante otro método.
Definamos la función escalón s(t):
=1
0)(ts
0
0
≥
<
t
t (4.66)
Cuando una excitación de tensión expresada por s(t) se aplica al puerto j de la
subred lineal ampliada, la respuesta del circuito definida como la corriente en el dominio del tiempo entrando en el puerto i se denotará como )(tsij . Consideremos
ahora una excitación de tensión genérica en banda base, g(t), aplicada al puerto j, que satisface:
)()( 1tgtg = ( )1tt ≤ (4.67)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
70
En el contexto del análisis MIHB, sólo son de interés los valores que g(t) toma en los instantes de tiempo de muestreo. Así g(t) puede aproximarse por:
[ ] [ ] [ ] .......)()()()()()()(1)()( 43332221 +−−−+−−−+−−≈ ttsttstgttsttstgttstgtg
[ ] )()()()(2
11 p
N
p
pp ttstgtgtg −−+= ∑=
− (4.68)
En el puerto i, la respuesta a la excitación de la subred lineal aumentada será una
corriente en banda base de la forma:
[ ] )()()()()0()(2
11 pij
N
p
ppAijLoij ttStgtgtgYtI −−+≈ ∑=
− (4.69)
donde )0(AijY es el valor correspondiente de la matriz de admitancias de la subred lineal
ampliada en DC. La muestra en el instante de muestreo n-ésimo ( Nn ≤≤1 ) de esa corriente vendrá dado por:
[ ] )()()()()0()(2
11 pnij
n
p
ppAijnLoij ttStgtgtgYtI −−+≈ ∑=
− (4.70)
donde el sumatorio es cero para n=1.
En el caso general, la excitación consiste en los )(tFo aplicados a los puertos de
polarización más los )(0 tU aplicados a los puertos del dispositivo. Asumiendo que el
índice j identifica a los puertos del dispositivo para Dnj ≤≤1 y a los puertos de
polarización para Bjb ≤≤ , )( nLoij tI queda:
[ ]∑ ∑∑= = =
− +−−+≈D Dn
j
n
j
pnij
n
p
pojpojojijnLoij ttStUtUtUYtI1 1 2
11 )()()()()0()(
[ ]∑ ∑∑= = =
− −−++Dn
j
B
j
pnij
n
p
pojpojojTij ttStFtFtFY1 1 2
11 )()()()()0(
(4.71)
Se asume que )(0 tF permanece constante a )( 10 tF para 1tt < . La expresión
(4.71) es la que sustituye a (4.65) para k=0. Formulación del sistema de resolución MIHB
Las corrientes lineal y no lineal en los puertos del dispositivo cumplen:
0)()( =+ titI L (4.72) O equivalentemente:
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
71
[ ]∑ =+ Ω
k
tj
kLkketWtI 0)()( (4.73)
Observemos que si las igualdades:
0)()()( =≅+ nknknLk tEtWtI (4.74)
donde los vectores )( nk tE son los errores complejos MIHB, se satisfacen para todo k y
n, y el espacio entre los instantes de muestro es suficientemente pequeño para que se cumpla el teorema de muestreo, entonces (4.74) es una condición suficiente para (4.73). De (4.74), el sistema de resolución queda:
,0)( =no tE
[ ] ,0)(Re =nk tE
[ ] 0)(Im =nk tE ( ).1;0: Nnk k ≤≤>Ω∀ (4.75)
Si denotamos P al número de productos de intermodulación de interés kΩ positivos,
(4.75) es un sistema de )12( += PNnN DT ecuaciones reales. 4.3.3 Solución del sistema MIHB Estrategia general
Debido a la peculiar naturaleza de los errores MIHB (X), una aproximación rápida a la solución de (4.75) puede encontrarse en seguida. Para un n dado y todo k que cumpla 0≥Ω k , vamos a colocar las partes real e imaginaria de los errores MIHB en un
vector de error real nE , y las partes real e imaginaria de los armónicos complejos en un
vector )( nk tX , al que nos referiremos como el vector de estado real del instante de
muestreo n-ésimo. (4.75) puede rescribirse de la siguiente forma:
[ ] 0,......,; 1 =−− Mnnnn SSSE )1( Nn ≤≤ (4.76)
Para un n dado, (4.76) puede verse como un sistema real de )12( +PnD
ecuaciones. Así, (4.76) puede resolverse como una secuencia de N sistemas independientes modificados de balance armónico, de tamaño )12( +PnD . Sin embargo, para Mn ≤≤1 , la información requerida para calcular las derivadas de la envolvente a través de (4.47) y (4.48) no está completamente disponible, con lo que las derivadas en esos instantes se calculan con un valor reducido de M ( 1' −= nM ). En particular, en 1t no hay información disponible sobre los estados del circuito en instantes de muestreo anteriores, ya que en esos instantes no se conocen las envolventes. Así, a 1t se lleva a cabo un análisis de balance armónico empezando desde cero armónicos, es decir, el circuito con fuentes puramente sinusoidales definidas por (4.32), con fasores constantes
)( 1tFk . Esto lleva a la determinación del primer vector de estado real, 1S . En los
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
72
siguientes pasos )1( >n , el sistema modificado de balance armónico se resuelve usando
el estado conocido 1−nS como punto de partida para la iteración.
Solución del sistema por métodos inexactos de Newton
La matriz jacobiana asociada al instante de muestreo n-ésimo se define por
n
nn
S
EJ
∂
∂=)( , y es una matriz cuadrada de dimensión )12( +PnD .
Según la discusión anterior, el algoritmo que vamos a ver se usará para calcular
)(nJ para todo n excepto 1=n . Para 1=n se lleva a cabo un análisis de balance armónico convencional, preferible al MIHB. Por lo tanto, en lo siguiente, sólo se considerará el caso 1>n .
Recordando que nE contiene las partes real e imaginaria de los errores MIHB
)( nk tE y que nS contiene las partes real e imaginaria de los armónicos )( nk tX , )(nJ
puede subdividirse en submatrices reales de la forma:
∂
∂
∂
∂∂
∂
∂
∂
=
))(Im(
))(Im(
))(Re(
))(Im())(Im(
))(Re(
))(Re(
))(Re(
)(,
ns
nk
ns
nk
ns
nk
ns
nk
n
sk
tX
tE
tX
tE
tX
tE
tX
tE
J (4.77)
donde la segunda fila se elimina para 0=k y la segunda columna se elimina para
0=s . Y sólo se consideran los valores de k y s tales que 0≥Ω k y 0≥Ω s . De la
expresión (4.74):
))((
)(
))((
)(
))((
)( ,
ns
nk
ns
nkL
ns
nk
tXG
tW
tXG
tI
tXG
tE
∂
∂+
∂
∂=
∂
∂ (4.78)
donde G denota la parte real o la imaginaria. Por otro lado, para 0≠k hacemos uso de (4.65) para obtener:
))((
)()()()(
))((
)(2
2
00,
ns
nk
k
ns
nkL
tXG
tU
d
Ydb
d
dYjaY
tXG
tI
kk∂
∂
Ω
Ω−
Ω
Ω−Ω=
∂
∂
Ω=ΩΩ=Ω
(4.79)
Y para 0=k y 1>n , obtenemos de (4.71):
∑= ∂
∂=
∂
∂ Dn
j ns
nj
ij
ns
niL
tXG
tUS
tXG
tI
1
00,
))((
)()0(
))((
)( (4.80)
Así, para evaluar la matriz jacobiana (4.77) a través de (4.78) sólo hay que
calcular las derivadas de )( nk tU y )( nk tW .
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
73
Como vimos, )( nk tU es el k-ésimo armónico de la función cuasiperiódica definida
por la primera ecuación de (4.58) y (4.59) a través de (4.52)-(4.55). Ahora vamos a calcular las derivadas parciales de ese armónico con respecto a las partes reales e imaginarias de )( ns tX . Primero introducimos las expansiones multidimensionales de
Fourier:
∑Ω
=∂
∂
p
tj
np
n
petCx
u)(,0
∑Ω
=∂
∂
p
tj
np
n
petCy
u)(,1
,1
∑Ω
=∂
∂
p
tj
np
n
petCy
u)(,2
,2
∑Ω
=∂
∂
p
tj
n
d
p
n
petCd
u)( (4.81)
Aplicándolo, tenemos:
[ ] [ ])()()()())(Re(
)(,1,1,0,0 nsknsksnsknsk
ns
nk tCtCjtCtCtX
tU+−+− −Ω++=
∂
∂
[ ] [ ]Tj
n
d
sk
Tj
n
d
sknsknsksss etCetCtCtC
Ω−
+
Ω−
−+− +++Ω− )()()()( ,2,22
[ ] [ ])()()()( ,2,20,1,10 nsknsknsknsk tCtCbtCtCa +−+− −+−+
[ ] [ ]Tj
n
d
sk
Tj
n
d
sknsknsksss etCetCTatCtCaj
Ω−
+
Ω−
−+− +−−Ω+ )()()()(2 0,2,20 (4.82)
[ ] [ ])()()()())(Im(
)(,1,1,0,0 nsknsksnsknsk
ns
nk tCtCtCtCjtX
tU+−+− +Ω−−=
∂
∂
[ ] [ ]Tj
n
d
sk
Tj
n
d
sknsknsksss etCetCjtCtCj
Ω−
+
Ω−
−+− −+−Ω− )()()()( ,2,22
[ ] [ ])()()()( ,2,20,1,10 nsknsknsknsk tCtCjbtCtCja +−+− −+−+
[ ] [ ]Tj
n
d
sk
Tj
n
d
sknsknsksss etCetCTjatCtCa
Ω−
+
Ω−
−+− −−+Ω− )()()()(2 0,2,20 (4.83)
En cuanto a las derivadas de )( nk tW que aparecen en (4.78), se calculan de
modo similar. Si las derivadas de las envolventes no aparecen en (4.42), (4.43) y (4.46), (4.76)
se reduciría a una secuencia de sistemas de balance armónico convencionales. Esta
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
74
situación sería equivalente a dejar 0== mm ba para todo m. Los cuatro primeros
términos en (4.82) y (4.83) son las expresiones que tendrían las derivadas si se hiciese un análisis convencional de balance armónico en un instante nt . El resto de términos
son las correcciones requeridas para explicar los acoplamientos entre los instantes de muestreo que son generados por la modulación de la envolvente.
Para resolver el sistema no lineal (4.76), una secuencia de iteraciones en el vector de estado incógnita se genera al ir resolviendo la ecuación de Newton en cada paso mediante un método de subespacios de Krylov, tales como GMRES. En este proceso, la gran parte de coste computacional usado se debe a la multiplicación de la matriz jacobiana por una secuencia de vectores reales, para construir los vectores base del subespacio de Krylov. Esto implica, debido a (4.78) y (4.79), que la operación básica a llevar a cabo es la evaluación de un gran número de sumatorios de la forma:
∑+
∂
∂+
∂
∂
s
s
ns
nk
s
ns
nk gtX
tUf
tX
tU
))(Im(
)(
))(Re(
)( (4.84)
(e igualmente con kW en vez de kU ), donde sf y sg son vectores reales de dimensión
Dn . El límite superior “+” indica que el sumatorio es extendido sólo a los vectores
armónicos s para los que 0≥Ω s . Observemos ahora que (4.82) y (4.83) pueden ser
reescritos de la forma sintetizada:
)()())(Re(
)(nsknsk
ns
nk tQtQtX
tU+− +=
∂
∂ (4.85)
[ ])()())(Im(
)(nsknsk
ns
nk tQtQjtX
tU+− +=
∂
∂ (4.86)
donde las expresiones completas de )( nsk tQ − y )( nsk tQ + pueden obtenerse fácilmente
por inspección de (4.82) y (4.83). Si introducimos las definiciones:
sss jgfz −− −= ( )0: <Ω∀ ss ,
,2 00 fz =
sss jgfz += ( )0: >Ω∀ ss , (4.87)
Entonces (4.84) toma la forma:
∑∑ −
+
=
∂
∂+
∂
∂
s
ssk
s
s
ns
nk
s
ns
nk zQgtX
tUf
tX
tU
))(Im(
)(
))(Re(
)( (4.88)
donde el último sumatorio está ahora extendido a todos los valores de s, y tiene la estructura de una convolución discreta, por lo que puede ser calculada eficientemente por la FFT.
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
75
La conclusión es que para cualquier instante de tiempo de muestreo dado, el sistema MIHB (4.76) puede ser resuelto por métodos de subespacios de Krylov con la misma eficiencia numérica que para análisis de balance armónico convencional.
4.4 Técnica alternativa
4.4.1 Introducción
Lo esencial del algoritmo propuesto es reducir drásticamente el tamaño de la matriz jacobiana necesaria en la solución del balance armónico, a través de la introducción de un nuevo concepto en la reducción de circuitos no lineales. Esto permite que el gran sistema de ecuaciones no lineales que representa la solución de balance armónico, se sustituya por un sistema mucho más pequeño de ecuaciones no lineales. Esta técnica fue propuesta por Emad Gal et al. [11]
La principal característica del sistema reducido, es que comparte con el original
las primeras derivadas con respecto al nivel de excitación RF en el punto de operación DC. Después se usa una técnica para trazar la solución de las ecuaciones reducidas y mapearla de vuelta al espacio de las ecuaciones del circuito original. Esto introduce una significativa ventaja computacional, porque el tamaño de la matriz jacobiana que necesita ser factorizada, es típicamente mucho menor que el del jacobiano original.
Recordemos brevemente la situación de partida. Por un lado, teníamos la ecuación de balance armónico escrita en el dominio de la frecuencia del siguiente modo:
0)()()( =+Ω++= VIVQjYVIVF S (4.89)
Y al aplicar el método de Newton-Raphson, vimos que el vector solución de los
armónicos se iba actualizando iterativamente del siguiente modo:
)()( )()(1)()1( jj
F
jj VFVJVV −+ −= (4.90)
Y si en vez de ese método, se usan métodos de continuación, se trata de aumentar el sistema en la ecuación X, con algún parámetro que podemos llamar α , para obtener un sistema auxiliar ),( αψ V . Aquí,α es elegido de modo que )0,(Vψ tiene solución trivial y )()1,( VFV =ψ . En nuestro caso, este sistema ampliado podemos
formularlo dividiendo el vector SI en dos vectores:
0)()(),( =+Ω+++≡ VIVQjYVIIV TDC ααψ (4.91)
donde TI es un vector que representa los armónicos debido a las fuentes RF en el
circuito, mientras que DCI representa las fuentes de polarización DC. Un método de
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
76
continuación típico consistiría en ir tomando pequeños incrementos de α y resolver el sistema ),( αψ V , hasta α =1. Aunque estos métodos evitan algunos problemas de convergencia, son computacionalmente caros. 4.4.2 Algoritmo propuesto El objetivo es reducir el tamaño del circuito original descrito por la ecuación (4.89). El conjunto de ecuaciones reducido se obtiene efectuando el cambio de variables
VUV ˆ→ en la ecuación (4.91), y luego premultiplicar por TU para obtener el nuevo sistema:
0)ˆ()ˆ(ˆˆˆˆ),(ˆ =+Ω+++≡ VUIUVUQjUVYIIV TT
TDC ααψ (4.92)
donde
YUUY T=ˆ , DC
T
DC IUI =ˆ , T
T
T IUI =ˆ (4.93)
En esta ecuación, xqHNU )12( +ℜ∈ es una base ortonormal para el subespacio expandido por los q primeros coeficientes en series de Taylor de V con respecto a α con q<< N(2H+1). En otras palabras, si la expansión en series de Taylor de V con respecto a α viene dado por:
∑∞
=
−=0
0 )(k
k
kAV αα (4.94)
entonces, colspan(U)=colspan( 0A , 1A ,..., 1−qA ) (4.95)
El sistema reducido obtenido con esta transformación mantiene los primeros q
coeficientes de la serie de Taylor del sistema original. Siendo más precisos, si V es expandido en serie de Taylor:
∑∞
=
−=0
0 )(ˆˆk
k
kAV αα (4.96)
entonces kk AUA ˆ= , k=0,...,q-1.
En este caso, aplicar el método de continuación al sistema en la ecuación (4.92) para trazar la trayectoria de la solución, requiere la factorización de la matriz jacobiana de dimensiones qxq:
UYdV
dQj
dV
dIUJ T
+Ω+=ˆ (4.97)
Esto supone una ventaja computacional significativa, ya que el tamaño de este jacobiano es, típicamente, mucho menor que el original. Las soluciones en el espacio de
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
77
las ecuaciones de balance armónico del circuito original se obtienen usando el mapeado
de la expresión VUV ˆ→ . 4.4.3 Cálculo de las derivadas Para calcular las derivadas, vamos a escribir (4.91) de la siguiente forma: 0)( =Φ+++ VYVII TDC α (4.98)
donde )()()( VIVQjV +Ω=Φ . Y ahora, expandimos V y )(VΦ en series de Taylor: ∑
=
=0
)(i
i
iAV αα , ∑=
=Φ0
))((i
i
iDV αα (4.99)
Claramente, en α =0, tenemos ))0((VΦ = 0D = )( 0AΦ . Sustituyendo en (4.98) y
estableciendo α =0 tenemos: 0)( 00 =Φ++ AYAI DC (4.100)
donde 0A vemos que es simplemente la solución DC. Para α =1:
011 =++ DYAIT (4.101) Ahora definimos )(αJ como:
dV
dJ
Φ=)(α , 00 )( == ααJJ (4.102)
Como dV
dA
d
dD
φ
α
φ
α
10
1 ===
, al sustituir obtenemos la siguiente relación para 1A :
TIAJY −=+ 10 )( (4.103)
Para calcular 2A ,..., nA , expandimos )(αJ como una serie de Taylor en α :
∑=
=0
)(k
k
kJJ αα (4.104)
y usamos esa expansión para escribir α
φ
d
d del siguiente modo:
αα
φ
d
dVJ
d
d= (4.105)
4. TÉCNICAS DE REDUCCIÓN DE CIRCUITOS
78
∑ ∑∑= =
−
=
− =1 1
1
0
1
i i
i
i
i
i
i
i
i iAJiD ααα (4.106)
Cogiendo la derivada n-ésima con respecto a α y haciendo α =0, obtenemos la siguiente relación recursiva:
∑−
=−−+=
1
10 )(
1 n
j
jnjnn AJjnn
AJD (4.107)
Así, tenemos:
∑−
=−−−=+
1
10 )(
1)(
n
j
jnjn AJjnn
AJY (4.108)
De las expresiones (4.103) y (4.108) queda claro que el cálculo de nAA ,...,1 sólo
requiere una descomposición LU de la matriz jacobiana, y cada derivada puede obtenerse con una sustitución hacia delante y hacia atrás. 4.4.4 Resumen del algoritmo Los principales pasos del algoritmo propuesto se resumen como sigue:
1) Empezando desde el punto de operación DC, se calculan las derivadas de los armónicos respecto a α , como se explica en el apartado de cálculo de derivadas.
2) Usando factorización QR se construye una base ortonormal
xqHNU )12( +ℜ∈ para las columnas de la matriz de derivadas de la expresión (4.95).
3) Con esa base, se formula un sistema reducido de ecuaciones no lineales,
como en las expresiones (4.92) y (4.93).
4) Se resuelven las ecuaciones (4.92) y (4.93) para V , incrementando poco a poco el valor de α ( 10 → ).
5) La solución del sistema original se obtiene en términos de V , usando
VUV ˆ→ , y se calcula el error residual. El orden de la reducción, q, se aumenta si el error residual es mayor que una tolerancia preespecificada.
Top Related