Initial commit

.idea/.gitignore

@@ -0,0 +1,8 @@
+# 디폴트 무시된 파일
+/shelf/
+/workspace.xml
+# 에디터 기반 HTTP 클라이언트 요청
+/httpRequests/
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml

.idea/SettlementPredictionCode_Step.iml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="inheritedJdk" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+</module>

.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

.idea/misc.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.9 (settle_prediction)" project-jdk-type="Python SDK" />
+</project>

.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/SettlementPredictionCode_Step.iml" filepath="$PROJECT_DIR$/.idea/SettlementPredictionCode_Step.iml" />
+    </modules>
+  </component>
+</project>

1_SP-11.csv

@@ -0,0 +1,80 @@
+Time,Settle,Surcharge
+0,0,1.5
+5,17.4,1.5
+7,23.9,1.5
+11,32.2,1.5
+14,41.7,1.5
+21,64.1,1.5
+28,72.5,1.5
+35,78.8,1.5
+42,93.3,1.5
+48,102.5,1.5
+53,108,3.002
+54,109.2,3.002
+55,110.4,3.002
+56,111.6,3.002
+59,117.3,3.002
+60,119.2,3.002
+61,121.1,3.002
+62,122.7,3.002
+67,130.2,3.002
+68,131.9,3.002
+69,133.6,3.002
+70,135.4,3.002
+74,141.4,3.002
+75,142.9,3.002
+76,144.4,3.002
+77,146.2,3.002
+80,149.2,3.002
+81,150.2,3.002
+82,151.2,3.002
+83,152.2,3.002
+91,162.8,3.002
+98,170,3.002
+105,177,3.002
+112,182.4,3.002
+115,185,3.002
+117,186.5,3.002
+118,187.3,3.002
+122,202.9,4.095
+124,210.5,4.095
+125,214.5,4.095
+126,218.6,4.095
+129,222.4,4.095
+130,223.7,4.095
+131,225,4.095
+132,226.3,4.095
+133,227.5,4.095
+136,231.7,4.095
+137,233.1,4.095
+138,234.5,4.095
+139,235.9,4.095
+140,237.3,4.095
+143,240.7,4.095
+147,245.5,4.095
+151,249.7,4.095
+154,252.8,4.095
+158,257.8,4.095
+161,261.1,4.095
+164,264.1,4.095
+168,268,4.095
+172,272.2,4.095
+175,275.5,4.095
+181,283.5,4.095
+192,293.5,4.095
+195,296.2,4.095
+199,301.3,4.095
+202,304.6,4.095
+209,311.1,4.095
+216,316,4.095
+223,322.3,4.095
+230,326.5,4.095
+237,331.6,4.095
+244,336.5,4.095
+251,341.2,4.095
+258,346.1,4.095
+266,350.9,4.095
+273,354,4.095
+280,356,4.095
+286,358,4.095
+294,360.9,4.095

main.py

@@ -0,0 +1,203 @@
+import numpy as np
+from scipy.optimize import least_squares
+
+# python library for visualization
+import matplotlib.pyplot as plt
+from matplotlib import rcParams
+
+# python library for data management
+import pandas as pd
+
+# generate a time-settlement curve for hyperbolic method
+def generate_data_hyper(px, pt):
+    return pt / (px[0] * pt + px[1])
+
+# error between regression and measurement
+def fun_hyper_linear(px, pt, py):
+    return pt / (px[0] * pt + px[1]) - py
+
+# Read .csv file using pandas
+data = pd.read_csv("1_SP-11.csv")
+
+# Set arrays for time and settlement
+time = data['Time'].to_numpy()
+settle = data['Settle'].to_numpy()
+surcharge = data['Surcharge'].to_numpy()
+
+# Set data range (in%) to use in the prediction
+start = 0
+end = 100
+
+# Find the data range (in data) to use in the prediction
+end_date = time[-1]
+pred_start_date = int(end_date * start / 100) # prediction start date
+pred_end_date = int(end_date * end / 100) # prediction end data
+
+# initialize the indices for start and end date
+start_index = -1
+end_index = -1
+
+# Find the index of the initial data
+count = 0
+for day in time: # time = [0, 1, 2, 3, ...]
+    count = count + 1
+    if day > pred_start_date:
+        start_index = count - 1
+        break
+
+# Find the index of the final data
+count = 0
+for day in time: # time = [... 100, 101, 104, ...]
+    count = count + 1
+    if day > pred_end_date:
+        end_index = count - 1
+        break
+
+# Set data for the prediction
+
+'''
+1단계 성토고 침하 예측
+'''
+
+#1단계에 해당하는 실측 데이터 범위 지정
+tm_1 = time[start:10]
+ym_1 = settle[start:10]
+
+# Set a list for the coefficient; here a and b
+x0 = np.ones(2)
+
+# declare a least square object
+res_lsq_hyper_linear_1 = least_squares(fun_hyper_linear, x0, args=(tm_1, ym_1))
+
+# Print the calculated coefficient
+print(res_lsq_hyper_linear_1.x)
+
+# Generate predicted settlement data
+settle_hyper_linear_1 = generate_data_hyper(res_lsq_hyper_linear_1.x, time)
+
+'''
+2단계 성토고 침하 예측
+'''
+
+# 2단계 실측 침하량 (2단계~최종)
+tm_2 = time[10:80]
+ym_2 = settle[10:80]
+# 2단계 실측 침하량 (2단계 구간만)
+tm_22 = time[10:37]
+ym_22 = settle[10:37]
+# 1단계 예측 침하량 (2단계 구간만)
+yp_2 = settle_hyper_linear_1[10:37]
+# 1단계 예측 침하량 (2단계 ~ 최종)
+yp_22 = settle_hyper_linear_1[10:80]
+
+# 2단계 보정 침하량 산정
+def fun_step_measured_correction(m, p):
+    return m - p
+
+step2_measured_correction = fun_step_measured_correction(ym_22, yp_2)
+
+# 2단계 t-ti 산정
+def fun_step_time_correction(t, ti):
+    return t - ti
+
+step2_time_correction = fun_step_time_correction(tm_2[0:69],tm_2[0])
+
+# 2단계 보정 침하량에 대한 예측 침하량 산정
+
+# declare a least square object
+# step2_time_correction[0:27]는 회귀분석 적용할 2단계 범위만 추출한 것
+res_lsq_hyper_linear_2 = least_squares(fun_hyper_linear, x0, args=(step2_time_correction[0:27], step2_measured_correction))
+
+# Print the calculated coefficient
+print(res_lsq_hyper_linear_2.x)
+
+# Generate predicted settlement data
+settle_hyper_linear_2 = generate_data_hyper(res_lsq_hyper_linear_2.x, step2_time_correction)
+
+# 2단계 침하곡선 작성
+def settlement_prediction_curve(m1, p1):
+    return m1 + p1
+
+step2_prediction_curve = settlement_prediction_curve(settle_hyper_linear_2, yp_22)
+
+# 2단계 보정 예측 침하량 산정
+def fun_step_prediction_correction(m2, p2):
+    return p2 + (m2[0] - p2[0])
+
+step2_prediction_correction = fun_step_prediction_correction(ym_2, step2_prediction_curve)
+
+'''
+3단계 성토고 침하 예측
+'''
+
+# 3단계 실측 침하량 (3단계~최종)
+tm_3 = time[37:end]
+ym_3 = settle[37:end]
+# 2단계 예측 침하량 (3단계 구간만)
+yp_3 = step2_prediction_correction[27:end]
+
+
+# 3단계 보정 침하량 산정
+step3_measured_correction = fun_step_measured_correction(ym_3, yp_3)
+
+# 3단계 t-ti 산정
+
+step3_time_correction = fun_step_time_correction(tm_3,tm_3[0])
+
+# 3단계 보정 침하량에 대한 예측 침하량 산정
+
+res_lsq_hyper_linear_3 = least_squares(fun_hyper_linear, x0, args=(step3_time_correction, step3_measured_correction))
+
+print(res_lsq_hyper_linear_3.x)
+
+settle_hyper_linear_3 = generate_data_hyper(res_lsq_hyper_linear_3.x, step3_time_correction)
+
+# 3단계 침하곡선 작성
+
+step3_prediction_curve = settlement_prediction_curve(settle_hyper_linear_3, yp_3)
+
+# 3단계 보정 예측 침하량 산정
+
+step3_prediction_correction = fun_step_prediction_correction(ym_3, step3_prediction_curve)
+
+'''
+그래프 작성
+'''
+
+# Set parameters for plotting
+rcParams['figure.figsize'] = (10, 10)
+
+# Subplot
+f, axes = plt.subplots(2,1)
+plt.subplots_adjust(hspace = 0.1)
+
+# draw surcharge data
+axes[0].plot(time, surcharge, color='black', label='surcharge height')
+
+axes[0].set_ylabel("Surcharge height (m)", fontsize = 17)
+axes[0].set_xlim(left = 0)
+
+# draw measured data
+axes[1].scatter(time, -settle, s = 50, facecolors='white', edgecolors='black', label = 'measured data')
+
+# draw predicted data
+axes[1].plot(time, -settle_hyper_linear_1, linestyle='--', color='red', label='original Hyperbolic_1')
+axes[1].plot(tm_2, -step2_prediction_correction, linestyle='--', color='blue', label='original Hyperbolic_2')
+axes[1].plot(tm_3, -step3_prediction_correction, linestyle='--', color='green', label='original Hyperbolic_3')
+
+# Set axes title
+axes[1].set_xlabel("Time (day)", fontsize = 17)
+axes[1].set_ylabel("Settlement (mm)", fontsize = 17)
+
+# Set min values of x and y axes
+axes[1].set_ylim(top = 0)
+axes[1].set_ylim(bottom = -1.5 * settle.max())
+axes[1].set_xlim(left = 0)
+
+# Set legend
+axes[1].legend(bbox_to_anchor = (0, 0, 1, 0), loc =4, ncol = 3, mode="expand",
+           borderaxespad = 0,  frameon = False, fontsize = 12)
+
+
+plt.savefig('1_SP_11_Rev.2.png', dpi=300)
+plt.show()