From 63d4de7076c650666b8bafb4952f5ae2f4fdd21d Mon Sep 17 00:00:00 2001 From: mayge Date: Mon, 22 Sep 2025 02:56:41 -0400 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=B7=BB=E5=8A=A0relax=E6=A8=A1?= =?UTF-8?q?=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core/relaxed_basisQR.py | 296 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 core/relaxed_basisQR.py diff --git a/core/relaxed_basisQR.py b/core/relaxed_basisQR.py new file mode 100644 index 0000000..b2591b8 --- /dev/null +++ b/core/relaxed_basisQR.py @@ -0,0 +1,296 @@ +import numpy as np +from core.sk_iter import generate_starting_poles +from scipy.linalg import block_diag +import skrf as rf +from skrf import VectorFitting +from core.freqency import auto_select + +class RelaxedBasicBasisQR: + def __init__(self,H,freqs,poles,weights=None,passivity=True,enforce_dc=True,fit_constant=True): + self.least_squares_rms_error = None + self.least_squares_condition = None + self.eigenval_condition = None + self.eigenval_rms_error = None + + self.dc_enforce = enforce_dc + self.fit_constant = fit_constant + + self.H = H + self.freqs = freqs + self.s = self.freqs * 2j * np.pi + self.P = len(poles) + self.poles = poles + self.Phi = self.generate_basis(self.s, self.poles) + self.A = self.matrix_A(self.poles) + self.B = self.vector_B(self.poles) + self.Cw,self.w0 = self.fit_denominator(self.H, weights=weights) + self.D = self.w0 + + self.Cr = None + + z = np.linalg.eigvals(self.A - self.B @ self.Cw) + p_next = -z + + if passivity: + self.next_poles = self.passivity_enforce(p_next) + else: + self.next_poles = p_next + # z = np.where(np.real(z) < 0, z, -np.conj(z)) # enforce LHP + # self.next_poles = np.sort_complex(z) + self.eigenval_condition = np.linalg.cond(self.A - self.B @ self.Cw) + self.eigenval_rms_error = np.sqrt(np.mean(np.abs(np.real(z) - np.real(poles))**2 + np.abs(np.imag(z) - np.imag(poles))**2)) + + self.Dt = self.eval_Dt_state_space() + self.delta = self.Dt / weights if weights is not None else self.Dt + pass + + def passivity_enforce(self,poles): + """enforce poles' real parts to be negative""" + enforced_poles = [] + for pole in poles: + if pole.real > 0: + pole = -np.conj(pole) + enforced_poles.append(pole) + return enforced_poles + + def eval_Dt_state_space(self): + """Return D(s_k)=C(s_k I - A)^(-1)B + D for all k (complex 1D array).""" + s = 1j * 2*np.pi * np.asarray(self.freqs, float).ravel() + A = np.asarray(self.A, np.complex128); n = A.shape[0] + B = np.asarray(self.B, np.complex128).reshape(n, 1) + C = np.asarray(self.Cw, float).reshape(1, n) + D = self.D + I = np.eye(n, dtype=np.complex128) + out = np.empty_like(s, dtype=np.complex128) + for k, sk in enumerate(s): + DS = D + (C @ np.linalg.inv(sk*I - A) @ B) + out[k] = DS[0, 0] + return out + + + def generate_basis(self,s, poles): + """Real basis of (15)-(16); returns Φ(s) and a layout for packing C.""" + cols = [] + i = 0 + while i < len(poles): + p = poles[i] + if p.real > 0: + raise ValueError("poles must be in the LHP") + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): + pc = poles[i+1] + phi1 = 1/(s - p) + 1/(s - pc) # eq (15)generate_basis + phi2 = 1j*(1/(s - p) - 1/(s - pc)) # eq (16) (fixed sign) + cols += [phi1, phi2] + i += 2 + else: + cols.append(1/(s - p)) + i += 1 + Phi = np.column_stack(cols).astype(np.complex128) + return Phi + + def matrix_A(self, poles): + def A_block(p): + if abs(p.imag) < 1e-14: + return np.array([[p.real]], float) # A_p = [ p ] + return np.array([[p.real, p.imag], # A_p = [[Re p, Im p], + [-p.imag, p.real]], float) # [-Im p, Re p]] + A = None; i = 0 + while i < len(poles): + p = poles[i] + Ab = A_block(p) + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): i += 2 + else: i += 1 + A = Ab if A is None else block_diag(A, Ab) + return A + + def vector_B(self, poles): + def B_block(p): + return np.array([[1.0]], float) if abs(p.imag)<1e-14 else np.array([[2.0],[0.0]], float) + B = None; i = 0 + while i < len(poles): + p = poles[i] + Bb = B_block(p) + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): i += 2 + else: i += 1 + B = Bb if B is None else np.vstack([B, Bb]) + return B + + def fit_denominator(self, H, weights=None): + """ + Solve formula (70) on the real basis Φ to obtain: + - d (real) → packs into C for this state's block structure + - gamma (complex) + Optional 'weights' (K,) apply row scaling: SK weighting if 1/|D_prev|. + """ + if weights is None: + weights = np.diag(np.ones(len(H), np.complex128)) + else: + weights = np.diag([1/res for res in weights]) + H = np.asarray(H, np.complex128).reshape(-1,1) + K, N = self.Phi.shape + one = np.ones((K, 1), np.complex128) + + psi = weights @ self.Phi + psi_w = np.hstack([weights@one, psi]) + Phi_w = np.hstack([one, self.Phi]) + + beta = np.linalg.norm((weights@H)) * len(H) + mean_row = (beta / K) * np.sum(Phi_w, axis=0) + + Hpsi_w = H * psi_w + + + + dc_tol = 1e-18 + has_dc = self.dc_enforce and self.freqs[0] < dc_tol + + if has_dc: + # Enforce DC response exactly: + k0 = int(np.argmin(np.abs(self.freqs))) + mask = np.ones(K, dtype=bool); mask[k0] = False + # Weighted rows for k≠0 + A_re = np.hstack([np.real(psi)[mask:], np.real(-Hpsi_w)]) + A_im = np.hstack([np.imag(psi)[mask:], np.imag(-Hpsi_w)]) + + else: + A_re = np.hstack([np.real(psi), np.real(-Hpsi_w)]) + A_im = np.hstack([np.imag(psi), np.imag(-Hpsi_w)]) + A_w_re = np.concatenate([np.zeros(N, float), np.real(mean_row).astype(float)]).reshape(1, -1) + + b_re = np.zeros_like(H) + b_im = np.zeros_like(H) + b_w_re = beta + + A = np.vstack([A_re, A_im, A_w_re]).astype(float) + Q,R = np.linalg.qr(A) + R22 = R[A_re.shape[1]//2:, A_re.shape[1]//2:] + Q2 = Q[:, A_re.shape[1]//2:] + # rown = np.linalg.norm(A, axis=1) + # rown = np.sqrt(rown) + # A = rown[:,None] * A + b = np.concatenate([b_re, b_im, [[b_w_re]]]).astype(float) + + x = np.linalg.inv(R22) @ (Q2.T @ b) + + self.least_squares_rms_error = np.sqrt(np.mean((Q2 @ R22 @ x - b)**2)) + self.least_squares_condition = np.linalg.cond(Q2 @ R22) + + Cw,w0 = self.vector_Cw(x,psi_w0=psi_w[:,0]) + return Cw,w0 + + def vector_Cw(self,x,psi_w0): + w0 = x[0] + C = x[1:] + return C.T,w0 + + def non_bias_Cr(self,w0): + A = np.asarray(self.Phi) + den = np.diag((w0 + self.Phi @ self.Cw.T).ravel()) + b = np.asarray(den) @ self.H.reshape(-1,1) + Cr, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) + return Cr + + def evaluate(self,freqs, w0): + s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() + phi = self.generate_basis(s, self.poles) + den = w0 + phi @ self.Cw.T + if self.Cr is None: + self.Cr = self.non_bias_Cr(w0=w0) + num = phi @ self.Cr + H = num / den + return H.ravel() + +if __name__ == "__main__": + start_point = 2 + network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") + K = 5 + H11,freqs = auto_select([network.y[i][0][0] for i in range(start_point,len(network.y))],network.f[start_point:],max_points=20) + poles = generate_starting_poles(2,beta_min=freqs[0]/1.1,beta_max=freqs[-1]*1.1) + + Dt_1 = np.ones((len(freqs),1),np.complex128) + # Levi step (no weighting): + basis = RelaxedBasicBasisQR(H11,freqs,poles=poles) + Dt = basis.Dt + poles = basis.next_poles + + print("Levi step (no weighting):") + print("A:",basis.A) + print("B:",basis.B) + print("C:",basis.Cw) + print("D:",basis.D) + print("next_pozles:",basis.next_poles) + print("Dt:",Dt, "norm:",np.linalg.norm(Dt)) + + # SK weighting (optional, after first pass): + least_squares_condition = [] + least_squares_rms_error = [] + eigenval_condition = [] + eigenval_rms_error = [] + for i in range(K): + basis = RelaxedBasicBasisQR(H11,freqs,poles=poles,weights=Dt) + Dt_1 = Dt + Dt = basis.Dt + poles = basis.next_poles + print(f"SK Iteration {i+1}/{K}") + print("A:",basis.A) + print("B:",basis.B) + print("C:",basis.Cw) + print("D:",basis.D) + print("z:",basis.next_poles) + print("Dt:",Dt) + print("Dt/Dt-1",np.linalg.norm(Dt) / np.linalg.norm(Dt_1)) + least_squares_condition.append(basis.least_squares_condition) + least_squares_rms_error.append(basis.least_squares_rms_error) + eigenval_condition.append(basis.eigenval_condition) + eigenval_rms_error.append(basis.eigenval_rms_error) + + # H11_evaluated = basis.evaluate_pole_residue(network.f[1:],poles,basis.C[0]) + H11_evaluated = basis.evaluate(network.f[start_point:], w0=basis.w0) + import matplotlib.pyplot as plt + fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) + + ax00 = axes[0][0] + ax00.plot(network.f[start_point:], np.abs([network.y[i][0][0] for i in range(start_point,len(network.y))]), 'o', ms=4, color='red', label='Samples') + ax00.plot(network.f[start_point:], np.abs(H11_evaluated), '-', lw=2, color='k', label='Fit') + ax00.plot(freqs, np.abs(H11), 'x', ms=4, color='blue', label='Input Samples') + ax00.set_title("Response i=0, j=0") + ax00.set_ylabel("Magnitude") + ax00.legend(loc="best") + + ax01 = axes[0][1] + ax01.set_title("Response i=0, j=0") + ax01.set_ylabel("Phase (deg)") + ax01.plot(network.f[start_point:], np.angle([network.y[i][0][0] for i in range(start_point,len(network.y))],deg=True), 'o', ms=4, color='red', label='Samples') + ax01.plot(network.f[start_point:], np.angle(H11_evaluated,deg=True), '-', lw=2, color='k', label='Fit') + ax01.plot(freqs, np.angle(H11,deg=True), 'x', ms=4, color='blue', label='Input Samples') + ax01.legend(loc="best") + + ax10 = axes[1][0] + ax10.plot(least_squares_condition, label='Least Squares Condition') + ax10.set_title("least_squares_condition") + ax10.set_ylabel("Magnitude") + ax10.legend(loc="best") + + ax11 = axes[1][1] + ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') + ax11.set_title("least_squares_rms_error") + ax11.set_ylabel("Magnitude") + ax11.legend(loc="best") + + ax20 = axes[2][0] + ax20.plot(eigenval_condition, label='Eigenvalue Condition') + ax20.set_title("eigenval_condition") + ax20.set_ylabel("Magnitude") + ax20.legend(loc="best") + + ax21 = axes[2][1] + ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') + ax21.set_title("eigenval_rms_error") + ax21.set_ylabel("Magnitude") + ax21.legend(loc="best") + fig.tight_layout() + plt.savefig(f"relaxed_basic_basis_QR.png") + + + +