search

convexengineering / gpfit

Fit posynomials to data
http://gpfit.readthedocs.io/en/latest/
MIT License
10 stars 7 forks source link

Empty array error #47

Open mayork opened 8 years ago

mayork commented 8 years ago

I've been trying to fit some compressor maps and am getting what appears to be a deterministic linear algebra error. I can get the error, then run the exact same code and it won't throw an error the second time. I'm not too familiar with gpfit so I'm not sure why this would occur.

Below is the error and some code which causes said error.

Traceback (most recent call last):
  File "/Users/mayork/Documents/GpGit/gpfit/gpfit/compressor_map_REAL_DATA_fitting.py", line 140, in <module>
    r,const = fit.fit(independent, dependent, 4, 'SMA')
  File "/Users/mayork/Documents/GpGit/gpfit/gpfit/fit.py", line 71, in fit
    bainit = max_affine_init(xdata, ydata, K)
  File "/Users/mayork/Documents/GpGit/gpfit/gpfit/max_affine_init.py", line 58, in max_affine_init
    if matrix_rank(X[inds, :]) < dimx + 1:
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 1543, in matrix_rank
    S = svd(M, compute_uv=False)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 1338, in svd
    _assertNoEmpty2d(a)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 222, in _assertNoEmpty2d
    raise LinAlgError("Arrays cannot be empty")
LinAlgError: Arrays cannot be empty
import numpy as np
import matplotlib.pyplot as plt
import fit

#variable to control if plotting occurs
PLOT = True

#-----------------------------------
#Compressor map
#make the data
#real data
##N = [.5, .6, .7, .75, .8, .85, .875, .9, .925, .95, .975, .985, 1, 1.025]
##pi= [[2.27690,2.21610,2.15520,2.00330,1.91210,1.66880],[3.21930,3.12810,3.00650,2.82410,2.61130,2.21610],[4.86100,4.58740,4.31380,4.16180,3.76660,3.18890],
##     [5.95550,5.65150,5.40830,5.19540,4.73940,4.00980],[7.68840,7.56680,7.32360,7.14120,6.62430,5.65150],[10.2117,9.72530,9.26930,8.93490,8.23560,7.20200],
##     [11.7318,11.3062,10.7590,10.2725,9.60370,8.38760],[13.6775,13.1303,12.5223,11.9446,11.2757,9.99890],[16.3529,15.8056,15.0760,14.4072,13.5255,12.1574],
##     [19.4235,18.9978,18.3594,17.4473,16.4745,14.9544],[23.4365,22.8284,22.0076,20.8523,19.8187,18.3898],[24.8350,24.3181,23.6493,22.5548,21.4604,20.1227],
##     [27.0239,26.2942,25.3822,24.1357,22.7372,21.5212],[29.0304,28.2704,27.2975,26.3246,24.5309,23.3453]]
##mbar = [[0.114208,0.118718,0.123241,0.126468,0.128404,0.129694],[0.155504,0.161956,0.165183,0.167763,0.170345,0.170345],[0.220674,0.229707,0.230998,0.232933,0.235515,0.236805],
##        [0.261969,0.274874,0.280681,0.282617,0.285843,0.287779],[0.345206,0.351659,0.364563,0.371015,0.376177,0.378759],[0.425216,0.436185,0.443283,0.448445,0.451026,0.452962],
##        [0.486515,0.496838,0.505227,0.509743,0.512324,0.512324],[0.560072,0.572332,0.580075,0.585237,0.586528,0.588463],[0.651697,0.667183,0.676216,0.682024,0.683314,0.684605],
##        [0.751710,0.770422,0.783972,0.792360,0.792360,0.792360],[0.871081,0.886569,0.896886,0.904634,0.907211,0.907862],[0.918187,0.934317,0.947862,0.958187,0.962057,0.963992],
##        [0.993033,0.997545,1.00142,1.00464,1.00593,1.00787],[1.05433,1.05820,1.06014,1.06337,1.06659,1.06659]]

##N = [.95, .985]
##pi= [[20.4125,20.4325,19.4235,18.9978,18.3594,17.4473,16.4745,14.9544],[25.8250,25.8350,24.8350,24.3181,23.6493,22.5548,21.4604,20.1227]]
##mbar = [[.552,.652,0.751710,0.770422,0.783972,0.792360,0.792360,0.79442360],[.7182,.8182,0.918187,0.934317,0.947862,0.958187,0.962057,0.964]]

uppi=[]
upm=[]
centerpi=26.19
centerm=1
for i in range(8):
    if i==0:
        uppi.extend([centerpi+.06])
        upm.extend([centerm-.217])
    if i==1:
        uppi.extend([centerpi+.063])
        upm.extend([centerm-.1368])
    if i==2:
        uppi.extend([centerpi+.042])
        upm.extend([centerm-.0618])
    if i==3:
        uppi.extend([centerpi])
        upm.extend([centerm])
    if i==4:
        uppi.extend([centerpi-.06])
        upm.extend([centerm+.0447])
    if i==5:
        uppi.extend([centerpi-.102])
        upm.extend([centerm+.0664])
    if i==6:
        uppi.extend([centerpi-.184])
        upm.extend([centerm+.0972])
    if i==7:
        uppi.extend([centerpi-.44])
        upm.extend([centerm+.0972])

uppi2=[]
upm2=[]
centerpi=17.9
centerm=.8
for i in range(8):
    if i==0:
        uppi2.extend([centerpi+.06])
        upm2.extend([centerm-.217])
    if i==1:
        uppi2.extend([centerpi+.063])
        upm2.extend([centerm-.1368])
    if i==2:
        uppi2.extend([centerpi+.042])
        upm2.extend([centerm-.0618])
    if i==3:
        uppi2.extend([centerpi])
        upm2.extend([centerm])
    if i==4:
        uppi2.extend([centerpi-.06])
        upm2.extend([centerm+.0447])
    if i==5:
        uppi2.extend([centerpi-.102])
        upm2.extend([centerm+.0664])
    if i==6:
        uppi2.extend([centerpi-.184])
        upm2.extend([centerm+.0972])
    if i==7:
        uppi2.extend([centerpi-.44])
        upm2.extend([centerm+.0972])

N=[1,.925]
pi=[uppi,uppi2]
mbar=[upm,upm2]
if PLOT == True:
#plot of data used in gpfit
    for i in range(len(N)):
        Nplot = N[i]*np.ones(len(mbar[0]))
        piplot = pi[i]
        mbarplot = mbar[i]
        plt.plot(mbarplot,piplot, '-r')
    plt.xlabel('Normalized Corrected Mass Flow')
    plt.ylabel('Fan Pressure Ratio')
    plt.title('E3 Fan Map')
    plt.show()

    for i in range(len(N)):
        Nplot = N[i]*np.ones(len(mbar[0]))
        piplot = pi[i]
        mbarplot = mbar[i]
        plt.plot(np.log(mbarplot),np.log(piplot), '-r')
    plt.xlabel('Log of Normalized Corrected Mass Flow')
    plt.ylabel('Log of Fan Pressure Ratio')
    plt.title('E3 Fan Map in Log Space')
    plt.show()

    for i in range(len(N)):
        Nplot = N[i]*np.ones(len(mbar[0]))
        invpiplot = np.ones(len(pi[i]))
        for j in range(len(pi[i])):
            invpiplot[j] = 1/(pi[i][j])
        mbarplot = mbar[i]
        plt.plot(np.log(mbarplot),np.log(invpiplot), '-r')
    plt.xlabel('Log of Normalized Corrected Mass Flow')
    plt.ylabel('Log of Fan Pressure Ratio')
    plt.title('E3 Fan Map in Log Space')
    plt.show()

#set up the data for the fit
Nfit = []
mbarfit = []
pifit = []

for i in range(len(N)):
    Nfit.extend(N[i]*np.ones(len(mbar[i])))
    for j in range(len(pi[i])):
        hold=pi[i][j]
        pifit.extend([1/hold])
    mbarfit.extend(np.divide(mbar[i],[1]))

#create the fit
independent = np.array([np.log(Nfit),np.log(mbarfit)])
dependent = np.log(pifit)
r,const = fit.fit(independent, dependent, 4, 'SMA')
print const

#plot the fit
nvec = np.linspace(.8, 1, 10)
mbarvec = np.linspace(.8,1,100)

for i in range(len(nvec)):
    N = nvec[i]
    pi=[]
    for j in range(len(mbarvec)):
        mbar = mbarvec[j]
        #fit to the tweaked data
        pi.extend([(0.282 * (N)**-3.56 * (mbar)**0.132
        + 9.75e-06 * (N)**-133 * (mbar)**49.7
        + 0.3 * (N)**-0.59 * (mbar)**-0.184
        + 0.306 * (N)**2.8 * (mbar)**0.0678)**(-1/.124)])
        #fit to the original data
##        pi.extend([(0.106 * (N)**-0.0299 * (mbar)**-0.129
##        + 0.119 * (N)**-0.0527 * (mbar)**-0.123
##        + 0.107 * (N)**0.028 * (mbar)**-0.146
##        + 0.102 * (N)**-0.0231 * (mbar)**-0.131
##        + 0.123 * (N)**-0.0233 * (mbar)**-0.131
##        + 0.136 * (N)**0.082 * (mbar)**-0.161)**(-1/.116)])

    plt.plot(mbarvec, pi, '-r')

#code for adding in the actual fan map data
##N = [.5, .6, .7, .75, .8, .85, .875, .9, .925, .95, .975, .985, 1, 1.025]
##pi= [[2.27690,2.21610,2.15520,2.00330,1.91210,1.66880],[3.21930,3.12810,3.00650,2.82410,2.61130,2.21610],[4.86100,4.58740,4.31380,4.16180,3.76660,3.18890],
##     [5.95550,5.65150,5.40830,5.19540,4.73940,4.00980],[7.68840,7.56680,7.32360,7.14120,6.62430,5.65150],[10.2117,9.72530,9.26930,8.93490,8.23560,7.20200],
##     [11.7318,11.3062,10.7590,10.2725,9.60370,8.38760],[13.6775,13.1303,12.5223,11.9446,11.2757,9.99890],[16.3529,15.8056,15.0760,14.4072,13.5255,12.1574],
##     [19.4235,18.9978,18.3594,17.4473,16.4745,14.9544],[23.4365,22.8284,22.0076,20.8523,19.8187,18.3898],[24.8350,24.3181,23.6493,22.5548,21.4604,20.1227],
##     [27.0239,26.2942,25.3822,24.1357,22.7372,21.5212],[29.0304,28.2704,27.2975,26.3246,24.5309,23.3453]]
##mbar = [[0.114208,0.118718,0.123241,0.126468,0.128404,0.129694],[0.155504,0.161956,0.165183,0.167763,0.170345,0.170345],[0.220674,0.229707,0.230998,0.232933,0.235515,0.236805],
##        [0.261969,0.274874,0.280681,0.282617,0.285843,0.287779],[0.345206,0.351659,0.364563,0.371015,0.376177,0.378759],[0.425216,0.436185,0.443283,0.448445,0.451026,0.452962],
##        [0.486515,0.496838,0.505227,0.509743,0.512324,0.512324],[0.560072,0.572332,0.580075,0.585237,0.586528,0.588463],[0.651697,0.667183,0.676216,0.682024,0.683314,0.684605],
##        [0.751710,0.770422,0.783972,0.792360,0.792360,0.792360],[0.871081,0.886569,0.896886,0.904634,0.907211,0.907862],[0.918187,0.934317,0.947862,0.958187,0.962057,0.963992],
##        [0.993033,0.997545,1.00142,1.00464,1.00593,1.00787],[1.05433,1.05820,1.06014,1.06337,1.06659,1.06659]]
##for i in range(len(N)):
##        Nplot = N[i]*np.ones(len(mbar[0]))
##        piplot = pi[i]
##        mbarplot = mbar[i]
##        plt.plot(mbarplot,piplot, '-b')

plt.xlabel('Normalized Corrected Mass Flow')
plt.ylabel('Pressure Ratio')
plt.title('Fan Map')
plt.show()
whoburg commented 8 years ago

Gpfit picks an initialization for the fitting process at random, so non determinism is expected.

whoburg commented 8 years ago

But, the empty array error is something we should track down.