Diferencia en resultados de Ecuaciones Lineales
Publicado por Sergio (6 intervenciones) el 09/12/2013 18:42:16
Que tal,
estoy resolviendo un sistema lineal Ax=B de matrices para encontrar solo 1columna de la matriz inversa A, esto lo realizo con varios métodos por ejemplo uso gesv y linsolve del paquete cvxopt, también calculo directamente la inversa de la matrix A para obtener la columna deseada. Sin embargo el resultado que obtengo es diferente dependiendo si la matriz la formo manualmente o si la formo de datos leidos de un archivo .cdf
Estos son los resultados que obtengo con la matriz creada manualmente:
Estos son los resultados que obtengo con la matriz creada con los datos leidos de un archivo .cdf:
si somos observadores notaremos que la parte real de dichas columnas es la única que presenta diferencia, esta diferencia es la que me genera problemas al ir avanzando en el problema que quiero resolver.
Fui descartando los posibles errores que pudiera tener y al parecer todo esta en orden, la matriz que se invierte (A) es la misma, la comprobación de que e realidad es la matriz inversa esta correcta (A*A^-1= matriz identidad), el número de condición de ambas matrices es igual (2428.4394189 para la matriz creada manualmente y 2428.51001029 para la matriz creada con el archivo .cdf). También comprobé la inversa de la matriz A en Matlab y resulta que es igual al primer caso.
les dejo el código y el archivo .cdf que leo y que estoy usando para tratar de sacar ese "fenómeno"...
Si alguien sabe el porque sucede esto se lo agradeceré.
Saludos y gracias
estoy resolviendo un sistema lineal Ax=B de matrices para encontrar solo 1columna de la matriz inversa A, esto lo realizo con varios métodos por ejemplo uso gesv y linsolve del paquete cvxopt, también calculo directamente la inversa de la matrix A para obtener la columna deseada. Sin embargo el resultado que obtengo es diferente dependiendo si la matriz la formo manualmente o si la formo de datos leidos de un archivo .cdf
Estos son los resultados que obtengo con la matriz creada manualmente:
1
2
3
4
5
6
7
[ 9.41e-03-j2.64e+00]
[ 1.33e-03-j2.72e+00]
[ 3.37e-03-j2.71e+00]
[ 3.37e-03-j2.71e+00]
[ 1.33e-03-j2.72e+00]
[ 9.41e-03-j2.67e+00]
[ 1.33e-03-j2.72e+00]
Estos son los resultados que obtengo con la matriz creada con los datos leidos de un archivo .cdf:
1
2
3
4
5
6
7
[ 5.04e-03-j2.64e+00]
[-3.03e-03-j2.72e+00]
[-1.00e-03-j2.71e+00]
[-1.00e-03-j2.71e+00]
[-3.03e-03-j2.72e+00]
[ 5.04e-03-j2.67e+00]
[-3.03e-03-j2.72e+00]
si somos observadores notaremos que la parte real de dichas columnas es la única que presenta diferencia, esta diferencia es la que me genera problemas al ir avanzando en el problema que quiero resolver.
Fui descartando los posibles errores que pudiera tener y al parecer todo esta en orden, la matriz que se invierte (A) es la misma, la comprobación de que e realidad es la matriz inversa esta correcta (A*A^-1= matriz identidad), el número de condición de ambas matrices es igual (2428.4394189 para la matriz creada manualmente y 2428.51001029 para la matriz creada con el archivo .cdf). También comprobé la inversa de la matriz A en Matlab y resulta que es igual al primer caso.
les dejo el código y el archivo .cdf que leo y que estoy usando para tratar de sacar ese "fenómeno"...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 03 16:06:37 2013
@author: Serge
"""
from numpy import array,concatenate,conjugate,real,imag,hstack,mean
from cvxopt import spmatrix,matrix,sparse,exp,spdiag,mul,div,normal
from cvxopt.umfpack import symbolic,numeric,solve,linsolve
from numpy.linalg import inv,solve,cond
from cvxopt.lapack import gesv,getrf
import numpy as np
#Ys=spmatrix([2+1j*-14.5,4-1j*8,6-1j*4,2-1j*2.5, 4+1j*8,5-1j*+17,3+1j*4,2+1j*5, 6+1j*4,3+1j*4,5+1j*-8.889, 2+1j*2.5,2+1j*5,2-1j*8.389],
#[0,1,2,3,0,1,2,3,0,1,2,0,1,3],[0,0,0,0,1,1,1,1,2,2,2,3,3,3])
ys1=[0-1j*30.5810, 0+1j*30.5810, 0-1j*27.0270, 0+1j*27.0270, 3.0425-1j*19.1702, -1.5212+1j*9.6314, -1.5212+1j*9.6314, 3.0425-1j*19.1702, -1.5212+1j*9.6314, -1.5212+1j*9.6314, 0+1j*27.0270, -1.5212+1j*9.6314, -1.5212+1j*9.6314, 3.0426-1j*73.2242, 0+1j*27.0270, 0+1j*30.5810, -1.5212+1j*9.6314, -1.5212+1j*9.6314, 3.0426-1j*49.7512, -0+1j*27.0270, 0-1j*27.0270]
ys2=[0,0,1,1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,6,6]
ys3=[0,5,1,4,2,4,5,3,4,5,1,2,3,4,6,0,2,3,5,4,6]
Ys=spmatrix(ys1,ys2,ys3,size=(7,7))
#print 'Ys:\n',Ys#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Vp=matrix((0.0+1j*0.0),(7,1))
Vp[0]=1.0+1j*0.0
#print Vp
gesv(matrix(Ys),Vp)
print Vp
B = matrix(Ys)
#B=np.matrix(Ys)
Zb=np.linalg.inv(B)
Zb=matrix(Zb)
iden=Ys*Zb
print '\n',iden
from numpy import array,concatenate,conjugate,real,imag,hstack
from math import sin,cos,pi,radians
from pylab import find
from cvxopt import spmatrix,matrix,sparse,exp,spdiag,mul,div
from cvxopt.umfpack import symbolic,numeric,solve
from cvxopt.lapack import gesvfrom time import time
from sys import argv
import numpy as np
Num_Bus=[];Tip_Bus=[];Vfi_Bus=[]; Afi_Bus=[];
CW_Bus=[];CVAR_Bus=[];GW_Bus=[];GVAR_Bus=[];
P_Rama=[];Q_Rama=[];R_Rama=[];X_Rama=[];
B_Rama=[];ITEMS=[];T_Rama=[];Tap_Tr=[];Bcomp=[];
tic=time()
#------------------Abriendo el archivo del caso de estudio
archivo=open('n1.cdf','r')
#------------------Leyendo la potencia base del sistema
Pot_base=archivo.readline()[30:37]
Pot_Base=float(Pot_base)
#print 'Potencia Base del Sistema:\n', '\t',Pot_Base
#------------------Leyendo el número de nodos del sistema
Numero_Nod=archivo.readline()
Numero_nod=Numero_Nod.find('ITEMS')
N_nod=Numero_Nod[Numero_nod-5:Numero_nod]
N_nod=int(N_nod)
#print 'Numero de nodos del Sistema:\n','\t',(N_nod)
#-------------------leyendo datos de buses
for i in xrange(N_nod):
dato=archivo.readline()
Num_Bus.append(int(dato[0:4]))#Numero del bus
Tip_Bus.append(int(dato[24:26]))#Tipo del bus
Vfi_Bus.append(float(dato[27:33]))#Voltaje final del bus
Afi_Bus.append(float(dato[33:40]))#Angulo final del bus
CW_Bus.append(float(dato[40:49]))#Carga en MW del bus
CVAR_Bus.append(float(dato[49:59]))#Carga en MVAR del bus
GW_Bus.append(float(dato[59:67]))#Generacion en MW del bus
GVAR_Bus.append(float(dato[67:75]))#Generacion en MVAR del bus
Bcomp.append(float(dato[114:123]))#Suceptancia de compensación
#---------------------Leyendo numero de ramas
dato2=archivo.readline()
N_ram=archivo.readline()
Numero_ram=N_ram.find('ITEMS')
N_Ram=N_ram[Numero_ram-5:Numero_ram]
N_Ram=int(N_Ram)
#print 'Numero de ramas del Sistema:\n','\t',N_Ram
#----------------------Formando arreglos de los datos leidos
Num_Bus=array(Num_Bus);Tip_Bus=array(Tip_Bus);Vfi_Bus=matrix(Vfi_Bus)
Afi_Bus=matrix(Afi_Bus);CW_Bus=array(CW_Bus);CVAR_Bus=array(CVAR_Bus)
GW_Bus=array(GW_Bus);GVAR_Bus=array(GVAR_Bus);Bcomp=array(Bcomp);
#Afi_Bus=Afi_Bus*pi/180
Afi_Bus[:]=0
#----------------------Obteniendo tipos de buses
#Bus tipo Carga
Tip_Carga=find(Tip_Bus==0)
Vfi_Bus[matrix(Tip_Carga)]=1.0
#Tip_Carga=array(Tip_Carga)
#print 'Numero de buses tipo carga',Tip_Carga+1
#Bus tipo PV
Tip_PV= find (Tip_Bus==2)
#print 'Numero de buses tipo PV',Tip_PV+1
#Bus tipo slack
Tip_Sl=find(Tip_Bus==3)
Tip_Sl=int(Tip_Sl)
#Buses con Geneadores
Bus_Gen=find(Tip_Bus!=0)
Num_Gen=len(Bus_Gen)
a=len(Bus_Gen)
NewGen=range(0,a)
Ndim=N_nod+len(Tip_Carga)
#------------------------Leyendo datos de ramas
for i in xrange(N_Ram):
dato3=archivo.readline()
P_Rama.append(int(dato3[0:4]))#Bus inicial
Q_Rama.append(int(dato3[5:9]))#Bus final
R_Rama.append(float(dato3[19:29]))#Resistencia de la rama
X_Rama.append(float(dato3[29:40]))#Reactancia de la rama
B_Rama.append(float(dato3[40:50]))#Susceptancia de la rama
T_Rama.append(int(dato3[18:19]))#Tipo de rama
Tap_Tr.append(float(dato3[76:82]))#Tap del transformador
#print B_Rama
archivo.close()
P_Rama=array(P_Rama); Q_Rama=array(Q_Rama);R_Rama=array(R_Rama);
X_Rama=array(X_Rama);B_Rama=array(B_Rama);T_Rama=array(T_Rama);
Tap_Tr=array(Tap_Tr);
#-------------------------Formando Ynodo
#Indices de lados p y q
pl=find(T_Rama!=1)
pt=find(T_Rama==1)
#Nodos p y q
PL=P_Rama[pl]
PT=P_Rama[pt]
QL=Q_Rama[pl]
QT=Q_Rama[pt]
#Indices de nodos con compensación
Ibc=find(Bcomp!=0)+1
#Susceptancia de compensación
Bcomp=1j*Bcomp[Ibc-1]
#Admitancia de la rama
Ylin=1/(R_Rama[pl]+X_Rama[pl]*1j)
Ytr=(1/(R_Rama[pt]+X_Rama[pt]*1j))/Tap_Tr[pt]
P1=concatenate((PL,PT))-1
P2=concatenate((QL,QT))-1
Y1=concatenate((-Ylin,-Ytr))
P3=concatenate((PL,QL,PT,QT,Ibc))-1
Y2=concatenate((Ylin+(1j*B_Rama[pl]/2),Ylin+(1j*B_Rama[pl]/2),\
Ytr/Tap_Tr[pt],Ytr*Tap_Tr[pt],Bcomp))
Ynodo=spmatrix(Y1,P1,P2,size=(N_nod,N_nod),tc='z')
Ynodo=Ynodo+Ynodo.T+spmatrix(Y2,P3,P3,size=(N_nod,N_nod),tc='z')
#print 'Ynodo:\n',Ynodo#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Zb1=np.linalg.inv(matrix(Ynodo))
Zb1=matrix(Zb1)
ide=Ynodo*Zb1
print ide
Vp1=matrix((0.0+1j*0.0),(7,1))
Vp1[0]=1.0+1j*0.0
#print Vp1
gesv(matrix(Ynodo),Vp1)
print Vp1
#print Ys-Ynodo
print 'Numero de condicion de Ys:',cond(matrix(Ys),'fro')
print 'Numero de condicion de Ynodo:',cond(matrix(Ynodo),'fro')
Si alguien sabe el porque sucede esto se lo agradeceré.
Saludos y gracias
- n1.zip(652,0 B)
Valora esta pregunta
0