I am trying to solve a DAE system using GEKKO. In the first part, I successfully solved the DAE. In the second part, I want to minimize one of my variables ('R'). Unfortunately, I cant find any solutions. Do you have any ideas what I did wrong? Here is my code:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
Tf = 150
T10 = 150.0
T20 = 150.0
T30 = 150.0
T40 = 150.0
T50 = 150.0
T60 = 150.0
Ta10 = Tf
Ta20 = Tf
Ta30 = Tf
Ta40 = Tf
Ta50 = Tf
Ta60 = Tf
T_f0 = Tf+1
Tw20 = 61
# Parameters:
# Stahl
m_steel = 7500 / 3600 # kg
cp = 500 # J/kg/K
# Luft
m_air = 100000 / 3600
cp_air = 10 # J/kg/K
alpha = 4500 # J/s
# Wasser
m_w = 0.15 # Strom durch HX
cp_w = 4200 # J/kg/K
# HX
UA = 120 # Wärmeübergangskoeffizient
Tw00 = 60 # Zulauftempratur
m = GEKKO() # create GEKKO model; "remote=False" for connection issues
R = m.MV(0.0, lb=0.0, ub=0.9) # Reflux
# m_w = m.MV(0.02, lb = 0.02, ub=1)
t = np.linspace(0, 30, 100)
m.time = t
# Charge
Q_c = np.zeros(100)
Q_c[2:40] = 300000
Q_charge = m.Param(value=Q_c)
# Discharge
Q_d = np.zeros(100)
Q_d[60:80] = -1
Q_discharge = m.Param(value=Q_d)
T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
T1.value = T10
T2.value = T20
T3.value = T30
T4.value = T40
T5.value = T50
T6.value = T60
Ta1.value = Ta10
Ta2.value = Ta20
Ta3.value = Ta30
Ta4.value = Ta40
Ta5.value = Ta50
Ta6.value = Ta60
T_f.value = T_f0
Tw1.value = Tw20
Tw2.value = Tw20
Tw12.value = Tw20
Q_Dis.value = 0
m.Equations([ # Storage
T1.dt() == (alpha * (T_f - T1)) / (m_steel * cp),
T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
T6.dt() == (alpha * (Ta5 - T6)) / (m_steel * cp),
# Calculate Temperature of air
Ta1.dt() == -(alpha * (Ta1 - T1)) / (m_air * cp_air),
Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
Ta6.dt() == -(alpha * (Ta6 - T6)) / (m_air * cp_air),
# HX
T_f == Q_Dis / (m_air * cp_air) + Ta6,
Tw2 == -Q_Dis / (m_w * cp_w) + Tw1,
# Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
Q_Dis == UA * (T_f - Tw1)/2 * Q_discharge,
# Nodes of the water system
Tw1 == ((m_w * (1 - (1-R))) * Tw12 + m_w * (1 - R) * Tw00) / ((m_w * (1 - (1 - R))) + m_w * (1 - R)),
Tw12 == Tw2
])
# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve()
# plot the prediction
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
plt.plot(m.time, T2.value, 'b-', linewidth=3)
plt.plot(m.time, T3.value, 'g-', linewidth=3)
plt.plot(m.time, T4.value, 'r--', linewidth=3)
plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)
plt.show()
# Optimize
m.options.IMODE = 6
# Tw2.UPPER = 250
# Tw2.LOWER = 110
R.Status = 1
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
T1.value = T1.value.value
T2.value = T2.value.value
T3.value = T3.value.value
T4.value = T4.value.value
T5.value = T5.value.value
T6.value = T6.value.value
Ta1.value = Ta1.value.value
Ta2.value = Ta2.value.value
Ta3.value = Ta3.value.value
Ta4.value = Ta4.value.value
Ta5.value = Ta5.value.value
Ta6.value = Ta6.value.value
T_f.value = T_f.value.value
Tw1.value = Tw1.value.value
Tw2.value = Tw2.value.value
Tw12.value = Tw12.value.value
Q_Dis.value = Q_Dis.value.value
m.Minimize(R)
m.solve(disp=True)
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
plt.plot(m.time, T2.value, 'b-', linewidth=3)
plt.plot(m.time, T3.value, 'g-', linewidth=3)
plt.plot(m.time, T4.value, 'r--', linewidth=3)
plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, R.value, 'r-', linewidth=3)
plt.show()
Kind Regars Cassian
I tried to simplifie my code as you can see here:
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
Tf = 150
T10 = 150.0
T20 = 150.0
T30 = 150.0
T40 = 150.0
T50 = 150.0
T60 = 150.0
Ta10 = Tf
Ta20 = Tf
Ta30 = Tf
Ta40 = Tf
Ta50 = Tf
Ta60 = Tf
T_f0 = Tf+1
Tw20 = 61
# Parameters:
# Stahl
m_steel = 7500 / 3600 # kg
cp = 500 # J/kg/K
# Luft
m_air = 100000 / 3600
cp_air = 10 # J/kg/K
alpha = 4500 # J/s
# Wasser
# R = 0.3 # Reflux
m_w = 0.15 # Strom durch HX
# m_w00 = 0.15 # Zustrom
# m_w22 = m_w * R # Abstrom
# m_w12 = m_w * (1 - R) # Reflux Strom
cp_w = 4200 # J/kg/K
# HX
UA = 120 # Wärmeübergangskoeffizient
Tw00 = 60 # Zulauftempratur
m = GEKKO() # create GEKKO model; "remote=False" for connection issues
R = m.MV(0.5, lb=0.49, ub=0.6)
# t = np.linspace(0, 30, 100)
t = np.linspace(0, 30, 50)
m.time = t
# Charge
# Q_c = np.zeros(100)
# Q_c[2:40] = 300000
Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)
# Discharge
# Q_d = np.zeros(100)
# Q_d[60:80] = -1
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)
# T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
T1, T6, Ta1, Ta6, T_f, Tw2, Tw1, Q_Dis = m.Array(m.Var, 8)
T1.value = T10
# T2.value = T20
# T3.value = T30
# T4.value = T40
# T5.value = T50
T6.value = T60
Ta1.value = Ta10
# Ta2.value = Ta20
# Ta3.value = Ta30
# Ta4.value = Ta40
# Ta5.value = Ta50
Ta6.value = Ta60
T_f.value = T_f0
Tw1.value = Tw20
Tw2.value = Tw20
# Tw12.value = Tw20
Q_Dis.value = 0
m.Equations([ # Storage
T_f / 1000 == (Q_Dis / (m_air * cp_air) + Ta6) / 1000,
Tw2 / 1000 == (-Q_Dis / (m_w * cp_w) + Tw1) / 1000,
# Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
Q_Dis / 10000 == UA * (T_f - Tw1) / 2 * Q_discharge / 10000,
# Nodes of the water system
Tw1 / 1000 == ((m_w * (R)) * Tw2 + m_w * (1 - R) * Tw00) / ((m_w * (R)) + m_w * (1 - R)) / 1000,
# Tw1 == R * Tw2 + R * Tw00 + Tw00,
# Tw12 / 1000 == Tw2 / 1000,
T1.dt() / 1000 == (alpha * (T_f - T1)) / (m_steel * cp) / 1000,
# T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
# T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
# T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
# T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
T6.dt() / 1000 == (alpha * (Ta1 - T6)) / (m_steel * cp) / 1000,
# Calculate Temperature of air
Ta1.dt() / 1000 == -(alpha * (Ta1 - T1)) / (m_air * cp_air) / 1000,
# Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
# Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
# Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
# Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
Ta6.dt() / 1000 == -(alpha * (Ta6 - T6)) / (m_air * cp_air) / 1000
])
# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve()
# plot the prediction
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)
plt.show()
# Optimize
m.options.IMODE = 6
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)
Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)
# Tw1.UPPER = 1250
# Tw2.LOWER = 10
R.Status = 1
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
T1.value = T1.value.value
# T2.value = T2.value.value
# T3.value = T3.value.value
# T4.value = T4.value.value
# T5.value = T5.value.value
T6.value = T6.value.value
Ta1.value = Ta1.value.value
# Ta2.value = Ta2.value.value
# Ta3.value = Ta3.value.value
# Ta4.value = Ta4.value.value
# Ta5.value = Ta5.value.value
Ta6.value = Ta6.value.value
T_f.value = T_f.value.value
Tw1.value = Tw1.value.value
Tw2.value = Tw2.value.value
# Tw12.value = Tw12.value.value
Q_Dis.value = 0 # Q_Dis.value.value
m.Minimize(R)
m.solve(disp=True)
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, R.value, 'r-', linewidth=3)
plt.show()
This code produces a solution, but not the one I expect or want haha. Some of my variables (Tw1 and Tw2) go against -∞ and the rest against +∞.
The simulation solves successfully, but the result is an unstable system.
The only objective is the minimize R
. Try adding another objective such as:
m.Minimize((Tw1-5)**2)
and also widen the bounds for R
such as R = m.MV(0.5, lb=0.2, ub=1.0)
. This gives a stable solution and an objective to guide the selection of R
.
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
Tf = 150
T10 = 150.0
T20 = 150.0
T30 = 150.0
T40 = 150.0
T50 = 150.0
T60 = 150.0
Ta10 = Tf
Ta20 = Tf
Ta30 = Tf
Ta40 = Tf
Ta50 = Tf
Ta60 = Tf
T_f0 = Tf+1
Tw20 = 61
# Parameters:
# Stahl
m_steel = 7500 / 3600 # kg
cp = 500 # J/kg/K
# Luft
m_air = 100000 / 3600
cp_air = 10 # J/kg/K
alpha = 4500 # J/s
# Wasser
# R = 0.3 # Reflux
m_w = 0.15 # Strom durch HX
# m_w00 = 0.15 # Zustrom
# m_w22 = m_w * R # Abstrom
# m_w12 = m_w * (1 - R) # Reflux Strom
cp_w = 4200 # J/kg/K
# HX
UA = 120 # Wärmeübergangskoeffizient
Tw00 = 60 # Zulauftempratur
m = GEKKO() # create GEKKO model; "remote=False" for connection issues
R = m.MV(0.5, lb=0.2, ub=1.0)
# t = np.linspace(0, 30, 100)
t = np.linspace(0, 30, 50)
m.time = t
# Charge
# Q_c = np.zeros(100)
# Q_c[2:40] = 300000
Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)
# Discharge
# Q_d = np.zeros(100)
# Q_d[60:80] = -1
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)
# T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
T1, T6, Ta1, Ta6, T_f, Tw2, Tw1, Q_Dis = m.Array(m.Var, 8)
T1.value = T10
# T2.value = T20
# T3.value = T30
# T4.value = T40
# T5.value = T50
T6.value = T60
Ta1.value = Ta10
# Ta2.value = Ta20
# Ta3.value = Ta30
# Ta4.value = Ta40
# Ta5.value = Ta50
Ta6.value = Ta60
T_f.value = T_f0
Tw1.value = Tw20
Tw2.value = Tw20
# Tw12.value = Tw20
Q_Dis.value = 0
m.Equations([ # Storage
T_f / 1000 == (Q_Dis / (m_air * cp_air) + Ta6) / 1000,
Tw2 / 1000 == (-Q_Dis / (m_w * cp_w) + Tw1) / 1000,
# Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
Q_Dis / 10000 == UA * (T_f - Tw1) / 2 * Q_discharge / 10000,
# Nodes of the water system
Tw1 / 1000 == ((m_w * (R)) * Tw2 + m_w * (1 - R) * Tw00) / ((m_w * (R)) + m_w * (1 - R)) / 1000,
# Tw1 == R * Tw2 + R * Tw00 + Tw00,
# Tw12 / 1000 == Tw2 / 1000,
T1.dt() / 1000 == (alpha * (T_f - T1)) / (m_steel * cp) / 1000,
# T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
# T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
# T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
# T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
T6.dt() / 1000 == (alpha * (Ta1 - T6)) / (m_steel * cp) / 1000,
# Calculate Temperature of air
Ta1.dt() / 1000 == -(alpha * (Ta1 - T1)) / (m_air * cp_air) / 1000,
# Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
# Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
# Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
# Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
Ta6.dt() / 1000 == -(alpha * (Ta6 - T6)) / (m_air * cp_air) / 1000
])
# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve()
# plot the prediction
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)
plt.show()
# Optimize
m.options.IMODE = 6
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)
Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)
# Tw1.UPPER = 1250
# Tw2.LOWER = 10
R.Status = 1
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
T1.value = T1.value.value
# T2.value = T2.value.value
# T3.value = T3.value.value
# T4.value = T4.value.value
# T5.value = T5.value.value
T6.value = T6.value.value
Ta1.value = Ta1.value.value
# Ta2.value = Ta2.value.value
# Ta3.value = Ta3.value.value
# Ta4.value = Ta4.value.value
# Ta5.value = Ta5.value.value
Ta6.value = Ta6.value.value
T_f.value = T_f.value.value
Tw1.value = Tw1.value.value
Tw2.value = Tw2.value.value
# Tw12.value = Tw12.value.value
Q_Dis.value = 0 # Q_Dis.value.value
m.Minimize(R)
m.Minimize((Tw1-5)**2)
m.solve(disp=True)
plt.figure(figsize=(8, 5))
plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
plt.subplot(3, 1, 3)
plt.plot(m.time, R.value, 'r-', linewidth=3)
plt.show()