diff --git a/optimizerO-series.py b/optimizerO-series.py
new file mode 100644
index 0000000..5182324
--- /dev/null
+++ b/optimizerO-series.py
@@ -0,0 +1,63 @@
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color_pool"]
+
+psgRange=np.arange(1e6,5e6,1e5)
+
+prediction = 2433827
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+temp2 = np.log(torShift)
+def tax(x):
+    return taxShift*(np.log(x)-temp1)/temp2
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+RNGk = 37.648854
+RNGb = 59397421.185785
+temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+def f1(x):
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+C2=1
+Cb2=1e8
+def f2(x):
+    return 1 - C2*x/Cb2 - 0.2
+
+C3=1
+Cb3=4e8
+def f3(x):
+    return 1 - C3*x/Cb3
+
+influenceFactor = np.array([0.21061,0.54848,0.24091])
+def f(x):
+    return f1(x)*influenceFactor[0] + f2(x)*influenceFactor[1] + f3(x)*influenceFactor[2]
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+print(result)
+
+for i in range(10):
+    C3=i*5
+    plt.plot(psgRange,f(psgRange),label='result, C3=%d'%C3,color=colorList[i])
+# plt.plot(psgRange,f1(psgRange),label='Ieco',color=colorList[1])
+# plt.plot(psgRange,f2(psgRange),label='Ienv',color=colorList[2])
+# plt.plot(psgRange,f3(psgRange),label='Isoc',color=colorList[3])
+    result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+    plt.plot(result.x,f(result.x),'o',color=colorList[i])
+    print(result.x,tax(result.x),f(result.x))
+
+plt.xlabel('total passengers')
+plt.ylabel('Optimized objective function')
+plt.legend(fontsize=8)
+plt.savefig('result/O-ser1.png',dpi=1024,)
+plt.show()
+
+print("Optimal total passengers:",result.x)
+print("Optimal taxation rate:",tax(result.x))
+print("Score",f(result.x))
\ No newline at end of file
diff --git a/stability-plot.py b/stability-plot.py
new file mode 100644
index 0000000..92bdeb6
--- /dev/null
+++ b/stability-plot.py
@@ -0,0 +1,24 @@
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+import csv
+
+with open('stability.csv', 'r') as f:
+    reader = csv.reader(f)
+    header = next(reader)
+    # next(reader)
+    data = [row for row in reader]
+    
+plt.figure(figsize=(6,5))
+plt.bar(range(len(data)), [np.log(np.abs(float(row[2]))) for row in data] , color=colorList[0])
+plt.plot([-1,len(data)],[0,0],color=colorList[2]) # x axis
+plt.margins(0.01)
+plt.xticks(range(0,len(data)), [row[3] for row in data], rotation=45)
+plt.xlabel('Parameters abbr.',labelpad=-3)
+plt.ylabel('Affect on Stability (log scale)')
+plt.title('Different Parameters\'s Affect on Stability')
+plt.savefig('result/stability-result.png',dpi=1024)
+plt.show()
\ No newline at end of file
diff --git a/stability.csv b/stability.csv
new file mode 100644
index 0000000..a9fcd89
--- /dev/null
+++ b/stability.csv
@@ -0,0 +1,12 @@
+name,algorithm,gradient
+Tax influences tourists,torShift**(1/taxShift),-38049184.08,TIT
+Passenger influence on revenue,the slope of the passenger-revenue relation curve,87730.21866,R-K
+Initial revenue,the interception of the passenger-revenue relation curve,-0.052009008,R-B
+Carbon footprint per person,C2,-230200.9189,Cpi
+Total Carbon footprint,Cb2,0.002691839,SCpi
+social effect per person,C3,-24062.72963,Spi
+maximum social effect bearage,Cb3,7.12E-05,SSpi
+iceberg contribution,iceberg,0.001602255,IC
+influence factor for ecomomy,influenceFactor[0],1574163.139,Weco
+influence factor for environment,influenceFactor[1],-432635.705,Wenv
+influence factor for society,influenceFactor[2],-100327.8061,Wsoc
diff --git a/stabilityO.py b/stabilityO.py
new file mode 100644
index 0000000..eed99a1
--- /dev/null
+++ b/stabilityO.py
@@ -0,0 +1,62 @@
+import numpy as np
+import json
+import matplotlib.pyplot as plt
+
+colorList = json.load(open('color/config.json','r'))["color"]
+
+psgRange=np.arange(1e6,5e6,1e5,dtype=np.float64)
+
+prediction = 2433827
+taxShift=0.03
+torShift=1-0.054
+curTaxationRate=1.0
+def tax(x):
+    temp1 = np.log(prediction/(torShift**(curTaxationRate/taxShift)))
+    temp2 = np.log(torShift)
+    return taxShift*(np.log(x)-temp1)/temp2
+def predictTotalPassengers(taxationRate):
+    return (prediction/(torShift**(curTaxationRate/taxShift)))*(torShift**(taxationRate/taxShift))
+RNGk = 37.648854
+RNGb = 59397421.185785
+def f1(x):
+    temp3=(curTaxationRate*(RNGk*(predictTotalPassengers(curTaxationRate))+RNGb))
+    return 5*((tax(x))*(RNGk*x+RNGb) / temp3 -1)+1
+
+C2=1
+Cb2=1e8
+iceberg = -0.2
+def f2(x):
+    return 1 - C2*x/Cb2 + iceberg
+
+C3=1
+Cb3=4e8
+def f3(x):
+    return 1 - C3*x/Cb3
+
+influenceFactor = np.array([0.21061,0.54848,0.24091],dtype=np.float64)
+def f(x):
+    return f1(x)*influenceFactor[0] + f2(x)*influenceFactor[1] + f3(x)*influenceFactor[2]
+
+
+from scipy.optimize import minimize_scalar
+
+result = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded')
+plt.plot(psgRange,f(psgRange),color=colorList[0],label='case 0')
+
+cur_x = iceberg
+cur_y = result.x
+tag = 1e-10
+alt_x = iceberg
+alt_y = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded').x
+while abs(alt_y-cur_y)<10:
+    iceberg += tag
+    alt_x = iceberg
+    alt_y = minimize_scalar(lambda x: -f(x),bounds=(np.min(psgRange),np.max(psgRange)),method='bounded').x
+    tag *= 2
+
+plt.plot(psgRange,f(psgRange),color=colorList[1],label='case 1')
+
+plt.legend()
+plt.show()
+print(cur_x,cur_y,alt_x,alt_y)
+print((alt_y-cur_y)/(alt_x-cur_x))
\ No newline at end of file