feat: 添加relax模式
This commit is contained in:
mayge
296
core/relaxed_basisQR.py
Normal file
296
core/relaxed_basisQR.py
Normal file
@@ -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")
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user