from scipy.optimize import leastsq import numpy as np import utils as ut # local copy from collections import OrderedDict # Copied from PT-4ZmeasVsModel script. # Can concatenate multiple measurement runs here. Make sure frequency and load arrays match. meas_mag = np.array([37.79])#, 54.01, 39.53, 19.94, 43.95, 48.36, 29.49, 30.12, 37.21, 36.31, 44.59, 29.47]) meas_pha = np.array([46.1])#, 28.9, 7.0, 18.3, 33.2, 20.8, 34.5, 9.9, 29.2, 20.5, 28.3, 20.9]) # Make complex nparray of impedances from mag & phase lists. This is the target # data to drive the model towards by adjusting PARAMS. target_data = meas_mag * np.exp(1j * np.radians(meas_pha)) # Fixed freq array for now meas_freqs = np.array([2.2e6] * len(meas_mag)) # Loads corresponding to frequencies and measurements # These are the "2.2 MHz Nominal Load Fixture" loads load_mag = np.array([44.6])#, 64.6, 81.5, 67.2, 56.7, 71.7, 56.1, 75.5, 59.6, 68.4, 61.8, 66.4]) load_pha = np.array([-48.4])#, -29.6, -50.6, -70, -39.8, -37.8, -58.6, -59.3, -48.6, -50.8, -39.9, -58.5]) zload = load_mag * np.exp(1j * np.radians(load_pha)) # Ordered dictionary of parameter names, initial guess values, and # min/max bounds. Use +/-np.inf for unbounded, set min = max to # fix a value. PINIT = 0 PBOUNDS = 1 BOUNDWEIGHT = 1e5 PARAMS = OrderedDict([ # name init_val (min, max) # ('Lpm', (700e-9, (10e-9, 1e-6))), # ('Rpm', (400e-3, (100e-3, 800e-3))), # ('Cpm', (40e-12, (20e-12, 100e-12))), ('phi', (0.0, (-20.0, 20.0))) ]) def model(w, params): L2 = 824e-9 # diff filter L C2 = 235e-12 # diff filter output C Lpm = 701e-9 # cable segment series L per meter Rpm = 1.00 # cable segment series R per meter Cpm = 46.5e-12 # cable segment shunt C per meter Ll2 = 280e-9 # CIC xfmr primary leakage L C3 = 6e-12 # CIC xfmr primary C Lm2 = 35.3e-6 # CIC xfmr magnetizing L N2z = 0.76563 # CIC xfmr pri:sec impedance ratio N2y = 1.0/N2z # same, admittance ratio Ll3 = 370e-9 # CIC xfmr secondary leakage L # Length of cable in meters LENGTH = 3.0 # Integer number of segments to split cable into SEGMENTS = 6 # Return the revised imaginary part of the immitance being adjusted add_xc = lambda C, im : im - 1.0/(w*C) add_xl = lambda L, im : im + w*L add_bc = lambda C, im : im + w*C add_bl = lambda L, im : im - 1.0/(w*L) # Get values to be optimized from the parameter list # Lps, Rps, Cps = [x * LENGTH / SEGMENTS for x in params[:3]] Lps = Lpm * LENGTH / SEGMENTS Rps = Rpm * LENGTH / SEGMENTS Cps = Cpm * LENGTH / SEGMENTS # Phase shift applied to final impedance # phi = params[3] phi = params[0] # Initialize z with load impedances, then add the effect of series & # shunt elements sequentially from load to measurement point (source) # (backward from the PT-4 calculation, which determines load from # measurement) z = zload # z = z + 0 # Uncomment this line to fix bug. z.imag = add_xl(Ll3, z.imag) # Bug in this line of code: # both z.imag AND zload.imag are # altered by this step, instead of # z.imag only y = 1.0 / z y *= N2y # CIC xfmr turns ratio y.imag = add_bl(Lm2, y.imag) # CIC xfmr magnetizing L y.imag = add_bc(C3, y.imag) # CIC xfmr primary C z = 1.0 / y z.imag = add_xl(Ll2, z.imag) # CIC xfmr primary leakage L y = 1.0/z # Iterate cable segments (input and output in y) for s in range(SEGMENTS): y.imag = add_bc(Cps, y.imag) z = 1.0 / y z.real += Rps z.imag = add_xl(Lps, z.imag) y = 1.0 / z y.imag = add_bc(C2, y.imag) # diff filter output C z = 1.0 / y z.imag = add_xl(L2, z.imag) # diff filter L # Phase shift the result z *= np.exp(1j * np.radians(phi)) # Return measured impedance return z def residuals(params, w, target_data): """ This is the error function minimized by leastsq. It should return the difference between the model and the data, not the square. :param params: tuple of variables to supply to model() :param w: radian frequency array :param target_data: impedance data to match :return: """ diff = target_data - model(w, params) # Split complex result into real & imag parts to minimize both zd = np.zeros(target_data.size*2, dtype=np.float64) zd[0:zd.size:2] = diff.real zd[1:zd.size:2] = diff.imag return zd def penalties(params, bounds, weight): # Build an array with number of parameters elements, 0 if the # parameter is within bounds and how far out if not: punishment = [np.fmin(x-lo, 0.0) + np.fmax(0.0, x-hi) for x, (lo, hi) in zip(params, bounds)] # Scale for how much it should hurt if you're out of bounds return weight * np.array(punishment) def b_residuals(params, w, target_data, bounds, weight): """ Create array of residuals including "penalty" elements when parameter min/max limits are exceeded. :param params: parameters adjusted to fit :param w: radian frequency array :param Z: target complex impedance array :return: array of residuals plus penalties """ ba = np.hstack((residuals(params, w, target_data), penalties(params, bounds, weight))) return ba def fit_model(target_data, f): # Convert f to angular freq w = 2 * np.pi * f # Extract the guess values out of the PARAMS values tuple init_val = [elem[PINIT] for elem in PARAMS.values()] # Extract bounds list of tuples from PARAMS values tuple bounds = [elem[PBOUNDS] for elem in PARAMS.values()] return leastsq(b_residuals, init_val, args=(w, target_data, bounds, BOUNDWEIGHT), full_output=1, maxfev=10000) # ======================================================================================= # Check input array sizes arrays_ok = len(target_data) == len(meas_freqs) and len(target_data) == len(zload) assert arrays_ok, "Data arrays not the same length" # Perform model optimization, extract results fitresult = fit_model(target_data, meas_freqs) fitparams = fitresult[0] infodict = fitresult[2] mesg = fitresult[3] ier = fitresult[4] # Print resulting model values for name, val in zip(PARAMS.keys(), fitparams): print('{0} = {1}'.format(name, ut.eng_notate(val))) # Print optimization info: print('Number of function calls: {0}'.format(str(infodict['nfev']))) if 1 <= ier <= 4: print('Found solution, ier = '.format(str(ier))) else: print('Solution not found. Message:') print(' {0}'.format(mesg))