Multi-Constellation GNSS Trilateration Code Error

Question

I’m trying to implement a multi-constellation trilateration code in Python to test an anti-spoofing mechanism. With reference of a single trilateration from this Github page, I attempted to execute a trilateration code capable of handling two constellations:

#Subset all satellites into groups of "num", and indicate how many constellations used
def satSubset(testdf, num):
    subsetLst = []
    for subset in combinations(testdf['Sat'], num):
        constLst = []
        for sat in subset:
            const = sat[0]
            constLst.append(const)
        constNum = len(set(constLst))
        subsetLst.append([subset, constNum])
    return subsetLst

#calculate satellite coordinates on transmission instead of reception
def referenceRotate(ECEF, indx, tau): 
    w = 7.2921151467e-5 #earth rotation in rad/sec
    xn = float(ECEF[indx][0])
    yn = float(ECEF[indx][1])
    zn = float(ECEF[indx][2])
    
    x_corrected = xn * np.cos(w*tau) + yn * np.sin(w*tau) 
    y_corrected = xn * -np.sin(w*tau) + yn * np.cos(w*tau)
    z_corrected = zn
    
    return x_corrected, y_corrected, z_corrected

#define function for optimization. This is the pseudorange/satellite coordinates equations matrix
def constResidual(sat_df, subset):
    def Fx_pos(guess):
        rows = []
        
        RINEX = sat_df['Sat'].tolist()
        PR_list = sat_df['PR'].tolist()
        clock_err_list = sat_df['CLK Bias'].tolist()
        ECEF_list = sat_df['ECEF'].tolist()
        
        sats_used = subset[0]
        constNum = subset[1]
        
        if constNum== 1:
            x, y, z, r = guess
            
            for sat in sats_used:
                indx = RINEX.index(sat)
                PR_corrected = float(PR_list[indx]) + float(clock_err_list[indx])*3e8
                tau = PR_list[indx]/3e8
            
                #print('SmPR of RINEX', RINEX[indx], 'is ', SmPR[indx])
                xn, yn, zn = referenceRotate(ECEF_list, indx, tau)
                #print('ECEF of RINEX', RINEX[indx], 'is ', [xn, yn, zn], '\n')

                val = np.sqrt((x - xn)**2 + (y - yn)**2 + (z - zn)**2) 
                eq = val - (PR_corrected - r)
                rows.append(eq)
            return rows
        
        elif constNum == 2:
            #For this project, 2 constellations will be used for multi-constellation purposes
            x, y, z, r_gps, r_gal = guess 
            
            for sat in sats_used:
                const = sat[0] #G or R
                indx = RINEX.index(sat)
                PR_corrected = float(PR_list[indx]) + float(clock_err_list[indx])*3e8
                tau = PR_list[indx]/3e8
            
                #print('SmPR of RINEX', RINEX[indx], 'is ', SmPR[indx])
                xn, yn, zn = referenceRotate(ECEF_list, indx, tau)
                #print('ECEF of RINEX', RINEX[indx], 'is ', [xn, yn, zn], '\n')
                val = np.sqrt((x - xn)**2 + (y - yn)**2 + (z - zn)**2) 
                
                if const == 'G':
                    eq = val - (PR_corrected - r_gps)
                    rows.append(eq)
                elif const == 'R':
                    eq = val - (PR_corrected - r_gal)
                    rows.append(eq)

            return rows
    return Fx_pos

def trilateration(sat_df, subset):
    #print('Satellite Subset: ', subset[0])
    initial_pos = list(np.zeros(subset[1] + 3))
    Fx_pos = constResidual(sat_df, subset)
    opt_pos = opt.least_squares(Fx_pos, initial_pos)
    
    return opt_pos, subset[0]

For satSubset, testdf is a dataframe consisting of available satellite RINEX ID, ECEF, and Pseudorange. For example:

For the most part, the residuals between the actual location and my calculations are acceptable. However, when there's only one satellite from one constellation, the least-squares optimization ignores it. For example, when using 5 GPS and 1 GLONASS, the locations derived are:

The two locations are the same despite using different GLONASS satellites. The trilateration code basically ignored constellations if there's only one satellite in it. When there's two or more, everything seems to be ok. What's wrong with my code?

Answer 1

The reason for my solution not working is because for multi-constellation trilateration to work, at least 2 satellites from each constellation is needed. In the above example, the GPS satellites will simply solve for x,y,z,r_gps, then solve for r_gal in the GLONASS satellite pseudorange equation.