/
Chronos
/
Control_lib
Обзор
Документация
Войти
/
Chronos
/
Control_lib
Код
Запросы
0
Задачи
Вики
Пакеты
0
Релизы
0
CI/CD
Аналитика
Безопасность
master
Operator_Norms.py
342 строки
11 KB
Chronos
Final version
03 май 2025, 14:49
03 май 2025, 14:49
0822167
Код
Авторство
О чём код?
import numpy as np import cvxpy as cp # # # Функция проверки выполнимости LMI def is_lmi_feasible_h_inf_continuous(A, B, C, gamma): n = A.shape[0] # Размер матрицы A m = B.shape[1] # Число входов P = cp.Variable((n, n), symmetric=True) # Переменная P LMI = cp.bmat([ [A.T @ P + P @ A + C.T @ C, P @ B], [B.T @ P, -gamma**2 * np.eye(m)] ]) constraints = [P >> np.eye(n) * 1e-6, LMI << -np.eye(n + m) * 1e-6] prob = cp.Problem(cp.Minimize(0), constraints) prob.solve(solver=cp.MOSEK, verbose=False, mosek_params={ 'MSK_DPAR_INTPNT_CO_TOL_PFEAS': 1e-7, 'MSK_DPAR_INTPNT_CO_TOL_DFEAS': 1e-7, 'MSK_DPAR_INTPNT_CO_TOL_REL_GAP': 1e-7 }) return prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] # Функция автоматического вычисления H∞-нормы def h_inf_norm_lmi_continuous(A, B, C, tol=1e-6): gamma_min = tol # Нижняя граница gamma_max = 1.0 # Начальная верхняя граница # Поиск верхней границы while not is_lmi_feasible_h_inf_continuous(A, B, C, gamma_max): gamma_max *= 2 # Бинарный поиск while gamma_max - gamma_min > tol: gamma = (gamma_min + gamma_max) / 2 if is_lmi_feasible_h_inf_continuous(A, B, C, gamma): gamma_max = gamma else: gamma_min = gamma return gamma # def is_lmi_feasible_h_inf_discrete(A, B, C, gamma): # n = A.shape[0] # Размер матрицы A # m = B.shape[1] # Число входов # P = cp.Variable((n, n), symmetric=True) # Переменная P # LMI = cp.bmat([ # [A.T @ P @A - P + C.T @ C, A.T@ P @ B], # [B.T @ P @A , -gamma**2 * np.eye(m)+B.T @ P @ B] # ]) # constraints = [P >> 0, LMI << 0] # prob = cp.Problem(cp.Minimize(0), constraints) # prob.solve(solver=cp.MOSEK, verbose=False) # return prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] # # def h_inf_norm_lmi_discrete(A, B, C, tol=1e-6, max_gamma=1e10): # gamma_min = tol # gamma_max = 1.0 # Initial upper bound # # # Find upper bound for gamma # while not is_lmi_feasible_h_inf_discrete(A, B, C, gamma_max): # gamma_max *= 2 # if gamma_max > max_gamma: # raise ValueError("H-infinity norm exceeds max_gamma. System may be unstable.") # # # Binary search to refine gamma # while gamma_max - gamma_min > tol: # gamma = (gamma_min + gamma_max) / 2 # if is_lmi_feasible_h_inf_discrete(A, B, C, gamma): # gamma_max = gamma # else: # gamma_min = gamma # # return gamma_max # def is_lmi_feasible_h_2(A, B, C, gamma): # n = A.shape[0] # # Переменная оптимизации: симметричная P > 0 # P = cp.Variable((n, n), symmetric=True) # # # Условия: P > 0, A P + P A^T + B B^T < 0, trace(C P C^T) <= gamma^2 # constraints = [ # P >> np.eye(n) * 1e-6, # A @ P + P @ A.T + B @ B.T << -np.eye(n) * 1e-6, # C @ P @ C.T - gamma**2 * np.eye(n) << 0 # ] # # # Задача: feasibility problem (нет целевой функции) # prob = cp.Problem(cp.Minimize(0), constraints) # prob.solve(solver=cp.SCS, verbose=False) # return prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] # # # def generalized_h_2_norm_lmi(A, B, C, tol=1e-6): # """ # Вычисляет обобщённую H2-норму системы через LMI. # x' = A x + B u, y = C x # # Параметры: # A, B, C: матрицы системы (numpy массивы). # gamma_min: начальная нижняя граница для бинарного поиска. # gamma_max: начальная верхняя граница для бинарного поиска. # tol: точность поиска. # # Возвращает: # Обобщённая H2-норма системы (число). # # Исключения: # ValueError: если система неустойчива. # """ # n = A.shape[0] # размерность состояния # m = B.shape[1] # число входов # p = C.shape[0] # число выходов # # # Проверка устойчивости системы # eigvals = np.linalg.eigvals(A) # if np.any(np.real(eigvals) >= 0): # raise ValueError("Система неустойчива, H2-норма не определена.") # # gamma_min = tol # Нижняя граница # gamma_max = 1.0 # Начальная верхняя граница # # # Поиск верхней границы # while not is_lmi_feasible_h_2(A, B, C, gamma_max): # gamma_max *= 2 # # # Бинарный поиск # while gamma_max - gamma_min > tol: # gamma = (gamma_min + gamma_max) / 2 # if is_lmi_feasible_h_2(A, B, C, gamma): # gamma_max = gamma # else: # gamma_min = gamma # return gamma # # # # # # # # # # # # def is_lmi_feasible_h_inf_control_continuous(A, Bu, Bv, C, D, gamma,return_sol=False): # n = A.shape[0] # Размерность состояния # m_v = Bv.shape[1] # Число возмущений (v) # m_u = Bu.shape[1] # Число управлений (u) # p = C.shape[0] # # Переменные оптимизации # Y = cp.Variable((n, n), symmetric=True) # P > 0 # Z = cp.Variable((m_u, n)) # K = Z * P^{-1} # # # # LMI из Bounded Real Lemma для синтеза # # LMI11 = A @ P + P @ A.T + Bu @ Z + Z.T @ Bu.T # # LMI12 = Bv # # LMI13 = (C @ P + D @ Z).T # # # # LMI21 = Bv.T # # LMI22 = -gamma * np.eye(m_v) # # LMI23 = np.zeros((m_v, p)) # # # # LMI31 = C @ P + D @ Z # # LMI32 = np.zeros((p, m_v)) # # LMI33 = -gamma * np.eye(p) # # # LMI = cp.bmat([ # # [A @ P + P @ A.T+ Bu @ Z + Z.T @ Bu.T+Bv@Bv.T,(C @ P + D @ Z).T], # # [C @ P + D @ Z,-gamma**2 * np.eye(p)] # # ]) # # LMI = cp.bmat([ # [A @ Y + Y @ A.T+ Bu @ Z + Z.T @ Bu.T, Bv, (C @ Y + D @ Z).T], # [Bv.T,-gamma * np.eye(m_v),np.zeros((m_v,p))], # [C @ Y + D @ Z,np.zeros((p,m_v)),-gamma * np.eye(p)] # ]) # # # Условия: P > 0 и LMI < 0 # constraints = [ # Y >> np.eye(n) * 1e-6, # LMI << -np.eye(n+m_v+p) * 1e-6 # ] # # # Решаем feasibility problem # prob = cp.Problem(cp.Minimize(0), constraints) # prob.solve(solver=cp.MOSEK, verbose=False, mosek_params={ # 'MSK_DPAR_INTPNT_CO_TOL_PFEAS': 1e-7, # 'MSK_DPAR_INTPNT_CO_TOL_DFEAS': 1e-7, # 'MSK_DPAR_INTPNT_CO_TOL_REL_GAP': 1e-7 # }) # if not return_sol: # return prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] # else: # return Y.value, Z.value # # def h_inf_norm_minimizing_control_lmi_continuous(A, Bu, Bv, C, D, tol=1e-6): # """ # Синтез Hinf-оптимального управления через LMI. # Система: # x' = A x + Bu u + Bv v # y = C x + D u # Управление: u = K x (K = Z * P^{-1}) # Возвращает (gamma_opt, K). # """ # n = A.shape[0] # Размерность состояния # m_v = Bv.shape[1] # Число возмущений (v) # m_u = Bu.shape[1] # Число управлений (u) # p = C.shape[0] # Число выходов (y) # # # gamma_min = tol # Нижняя граница # gamma_max = 1.0 # Начальная верхняя граница # # # Поиск верхней границы # while not is_lmi_feasible_h_inf_control_continuous(A, Bu, Bv, C, D, gamma_max): # gamma_max *= 2 # # while gamma_max - gamma_min > tol: # gamma = (gamma_min + gamma_max) / 2 # if is_lmi_feasible_h_inf_control_continuous(A, Bu, Bv, C, D, gamma): # gamma_max = gamma # else: # gamma_min = gamma # # P_opt,Z_opt = is_lmi_feasible_h_inf_control_continuous(A, Bu, Bv, C, D, gamma,return_sol=True) # # Вычисляем матрицу обратной связи K = Z * P^{-1} # K = Z_opt @ np.linalg.inv(P_opt) # return gamma, K # # # # # # # # # # # # # # # # # # # # # # # # # def is_lmi_feasible_h_2_control(A, Bu, Bv, C, D, gamma,return_sol=False): # n = A.shape[0] # Размерность состояния # m_v = Bv.shape[1] # Число возмущений (v) # m_u = Bu.shape[1] # Число управлений (u) # p = C.shape[0] # # Переменные оптимизации # # Переменные оптимизации # P = cp.Variable((n, n), symmetric=True) # P > 0 (n×n) # Z = cp.Variable((m_u, n)) # K = Z * P^{-1} (m_u×n) # # # Условия LMI для обобщённой H₂-нормы # # Условие 1: P > 0 # constraint_p = P >> 1e-6 * np.eye(n) # # # Условие 2: [P (C P + D Z)ᵀ # # (C P + D Z) γ² I] > 0 # LMI_peak = cp.bmat([ # [P, (C @ P + D @ Z).T], # [C @ P + D @ Z, gamma**2 * np.eye(p)] # ]) # constraint_peak = LMI_peak >> 1e-6 * np.eye(n + p) # Размер: (n+p)×(n+p) # # # Условие 3: A P + P Aᵀ + Bv Bvᵀ + Bu Z + Zᵀ Buᵀ < 0 # LMI_stability = ( # A @ P + P @ A.T # + Bv @ Bv.T # + Bu @ Z # + Z.T @ Bu.T # ) # constraint_stability = LMI_stability << -1e-6 * np.eye(n) # Размер: n×n # # # Собираем все условия # constraints = [ # constraint_p, # constraint_peak, # constraint_stability # ] # # # Решаем feasibility problem # prob = cp.Problem(cp.Minimize(0), constraints) # prob.solve(solver=cp.SCS, verbose=False) # if not return_sol: # return prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE] # else: # return P.value, Z.value # # def generalized_h_2_norm_minimizing_control_lmi(A, Bu, Bv, C, D, tol=1e-6): # """ # Синтез обобщённого H2-оптимального управления через LMI. # Система: # x' = A x + Bu u + Bv v # y = C x + D u # Управление: u = K x (K = Z * P^{-1}) # Возвращает (gamma_opt, K). # """ # n = A.shape[0] # Размерность состояния # m_v = Bv.shape[1] # Число возмущений (v) # m_u = Bu.shape[1] # Число управлений (u) # p = C.shape[0] # Число выходов (y) # # gamma_min = tol # Нижняя граница # gamma_max = 1.0 # Начальная верхняя граница # # # Поиск верхней границы # while not is_lmi_feasible_h_2_control(A, Bu, Bv, C, D, gamma_max): # gamma_max *= 2 # # # # # while gamma_max - gamma_min > tol: # gamma = (gamma_min + gamma_max) / 2 # if is_lmi_feasible_h_2_control(A, Bu, Bv, C, D, gamma): # gamma_max = gamma # else: # gamma_min = gamma # # P_opt,Z_opt = is_lmi_feasible_h_2_control(A, Bu, Bv, C, D, gamma,return_sol=True) # # # Вычисляем матрицу обратной связи K = Z * P^{-1} # K = Z_opt @ np.linalg.inv(P_opt) # return gamma, K