import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
# 1. 단계: 하중을 입력으로 받는 범용 PINN 모델 (beam의 처짐 예측)
class PINN_MultiLoad(nn.Module):
def __init__(self):
super(PINN_MultiLoad, self).__init__()
self.net = nn.Sequential(
nn.Linear(2, 64), nn.Tanh(), # 입력: (x, q) 2차원. 위치와 하중값
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, 1)
)
def forward(self, xq): # xq: (batch, 2) - [x, q]
return self.net(xq)
# 2. 단계: 물리 법칙 및 경계 조건 포함 손실 함수. Euler-Bernoulli Beam Theory 고려.
def compute_loss(model, xq_train, w_train, xq_physics, E, I, L, w_normalize, device):
# (1) Data Loss: 계측 데이터와의 오차
w_pred = model(xq_train)
loss_data = torch.mean((w_pred - w_train)**2) # 처짐 예측과 실제 처짐 간의 MSE
# (2) Physics Loss: PDE 검증
xq_physics_grad = xq_physics.clone().requires_grad_(True)
w = model(xq_physics_grad)
# 자동 미분 (x에 대해서만). Euler-Bernoulli 보 이론에 따르면, d⁴w/dx⁴ = q/(E*I) 이므로, 4차 미분이 필요함. 그래야, 모델이 다양한 하중 조건에서 물리 법칙인 d⁴w/dx⁴ = q/(E*I)를 만족하도록 학습할 수 있음.
dw_dx = torch.autograd.grad(w, xq_physics_grad, torch.ones_like(w), create_graph=True)[0][:, 0:1] # x에 대한 1차 미분. 처짐의 변화율. 기울기
d2w_dx2 = torch.autograd.grad(dw_dx, xq_physics_grad, torch.ones_like(dw_dx), create_graph=True)[0][:, 0:1] # x에 대한 2차 미분. 곡률. 휨
d3w_dx3 = torch.autograd.grad(d2w_dx2, xq_physics_grad, torch.ones_like(d2w_dx2), create_graph=True)[0][:, 0:1] # x에 대한 3차 미분. 곡률의 변화율. 휨의 변화
d4w_dx4 = torch.autograd.grad(d3w_dx3, xq_physics_grad, torch.ones_like(d3w_dx3), create_graph=True)[0][:, 0:1] # x에 대한 4차 미분. 휨의 변화율. 보의 강성에 대한 저항력
# q는 입력에서 추출 (역정규화)
q_batch = xq_physics[:, 1:2] * 20000.0 - 10000.0 # 정규화. q를 실제 하중 범위로 변환. q는 음수이므로 -10000에서 시작하여 20000 범위로 확장됨. 예: q=0 → -10000 N/m, q=1 → 10000 N/m. 만약, q가 -15000에서 -5000 범위라면, q_batch = xq_physics[:, 1:2] * 10000.0 - 15000.0로 조정해야 함.
# 정규화된 PDE: d⁴w̄/dx̄⁴ = (q*L⁴)/(E*I*w_max)
q_normalized = (q_batch * L**4) / (E * I * w_normalize) # q를 정규화하여 PDE의 잔차 계산. 참고로 PDE 약어는 Partial Differential Equation의 약자로 부분 미분 방정식을 의미함. 보의 처짐과 하중 사이의 관계를 나타내는 PDE의 잔차를 계산하여 모델이 물리 법칙을 얼마나 잘 만족하는지 평가함.
pde_residual = d4w_dx4 - q_normalized # 모델이 예측한 분포하중 q에 대한 4차 미분과 실제 q에 대한 정규화된 PDE의 차이
loss_physics = torch.mean(pde_residual**2)
# (3) Boundary Condition Loss: 다양한 q에 대해 경계 조건
# 샘플링: 여러 q 값에서 경계 조건 테스트
num_bc_samples = 10
q_bc = torch.linspace(0, 1, num_bc_samples).view(-1, 1).to(device)
# x=0에서 w=0. 보의 양끝은 고정되어 있다고 가정
xq_bc_0 = torch.cat([torch.zeros(num_bc_samples, 1, device=device), q_bc], dim=1)
w_bc_0 = model(xq_bc_0)
# x=1에서 w=0. 보의 양끝은 고정되어 있다고 가정
xq_bc_L = torch.cat([torch.ones(num_bc_samples, 1, device=device), q_bc], dim=1)
w_bc_L = model(xq_bc_L)
loss_bc = torch.mean(w_bc_0**2) + torch.mean(w_bc_L**2) # 양 끝단에서 처짐이 0이 되도록 강제
# 전체 손실
weight_data = 1.0 # 데이터 손실의 중요성
weight_physics = 3.0 # 물리PDE의 중요성
weight_bc = 50.0 # 경계 조건의 중요성
return weight_data * loss_data + weight_physics * loss_physics + weight_bc * loss_bc
# 3. 단계: 다양한 하중 조건에 대한 데이터 생성
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(f"사용 장치: {device}")
L, E, I = 10.0, 210e9, 0.0001
# 여러 하중 조건 생성 (q: -15000 ~ -5000 N/m 범위)
q_values = np.linspace(-15000, -5000, 10) # 10가지 하중 조건
x_positions = np.linspace(0, L, 21) # 각 하중마다 21개 위치
# 전체 데이터 생성
xq_train_list = []
w_train_list = []
for q in q_values:
for x in x_positions[::2]: # 10개 포인트만 사용
x_norm = x / L
q_norm = (q + 10000) / 20000 # q를 [0,1] 범위로 정규화
# 해석해
w_true = (q * x / (24 * E * I)) * (L**3 - 2 * L * x**2 + x**3)
xq_train_list.append([x_norm, q_norm])
w_train_list.append([w_true])
xq_train_np = np.array(xq_train_list)
w_train_np = np.array(w_train_list)
# 정규화
w_normalize = np.abs(w_train_np).max()
w_train_norm = w_train_np / w_normalize
print(f"데이터 정규화: x_scale={L}, w_scale={w_normalize:.6e}, q_range=[-15000, -5000]")
print(f"훈련 데이터 수: {len(xq_train_list)}")
# 텐서 변환
xq_train = torch.FloatTensor(xq_train_np).to(device)
w_train = torch.FloatTensor(w_train_norm).to(device)
# 물리 법칙 검증용 포인트 (여러 q 값)
x_phys = torch.linspace(0, 1, 50).view(-1, 1)
q_phys = torch.rand(50, 1) # 랜덤 q 값
xq_physics = torch.cat([x_phys, q_phys], dim=1).to(device)
# 4. 단계: 모델 및 학습 설정
model = PINN_MultiLoad().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=5e-4)
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=500)
# 5. 단계: 학습 루프
print("범용 하중 조건 PINN 학습 시작...")
epochs = 20000
for epoch in range(epochs + 1):
optimizer.zero_grad()
loss = compute_loss(model, xq_train, w_train, xq_physics, E, I, L, w_normalize, device)
loss.backward()
optimizer.step()
scheduler.step(loss)
if epoch % 1000 == 0:
current_lr = optimizer.param_groups[0]['lr']
print(f"Epoch {epoch}: Total Loss {loss.item():.6f}, LR={current_lr:.6e}")
# 6. 단계: 테스트 - 여러 하중 조건에서 예측
model.eval()
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
test_loads = [-15000, -10000, -7500, -5000]
for idx, (ax, q_test) in enumerate(zip(axes.flat, test_loads)):
with torch.no_grad():
# 테스트 입력 생성
x_test = torch.linspace(0, 1, 100).view(-1, 1)
q_test_norm = (q_test + 10000) / 20000
q_test_tensor = torch.full((100, 1), q_test_norm)
xq_test = torch.cat([x_test, q_test_tensor], dim=1).to(device)
# 예측
w_pred_norm = model(xq_test).cpu().numpy()
w_pred = w_pred_norm * w_normalize
# 해석해
x_full = np.linspace(0, L, 100)
w_exact = (q_test * x_full / (24 * E * I)) * (L**3 - 2 * L * x_full**2 + x_full**3)
ax.plot(x_full, w_exact, 'r-', linewidth=2, label='Exact Solution')
ax.plot(x_full, w_pred, 'b--', linewidth=2, label='PINN Prediction')
ax.set_xlabel('Position x (m)')
ax.set_ylabel('Deflection w (m)')
ax.set_title(f'Load q = {q_test} N/m')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.savefig('pinn_multiload_results.png', dpi=150)
plt.show()
# 7. 모델 사용 예시
print("\n모델 사용 예시")
def predict_deflection(model, x_meter, q_load, L, w_normalize, device):
"""임의의 위치와 하중에서 처짐 예측"""
x_norm = x_meter / L
q_norm = (q_load + 10000) / 20000
xq_input = torch.tensor([[x_norm, q_norm]], device=device)
with torch.no_grad():
w_norm = model(xq_input).cpu().item()
w_real = w_norm * w_normalize
return w_real
# 테스트
test_cases = [
(5.0, -10000), # 중앙, 원래 하중
(5.0, -12000), # 중앙, 더 큰 하중
(3.0, -8000), # 왼쪽, 작은 하중
]
for x, q in test_cases:
w_pred = predict_deflection(model, x, q, L, w_normalize, device)
print(f"위치 {x}m, 하중 {q} N/m → 예측 처짐: {w_pred:.6f} m")
