#!/usr/bin/env python
# -*- coding: utf-8 -*-
#

from __future__ import division
from __future__ import print_function

import matplotlib.pyplot as plt
from pyomo.core.kernel.expr import exp
from pyomo.core.kernel.numvalue import value
from pyomo.dae import *
#: pyomo imports
from pyomo.environ import *
from pyomo.opt import SolverFactory, SolverStatus

__author__ = "<NAME>"  #: May 2018

def main():
    with_plots = False
    #: Number of finite elements
    nfe = 100
    #: Number of collocation points
    ncp = 3
    m = ConcreteModel()

    m.nfe = nfe
    m.ncp = ncp

    m.t = ContinuousSet(bounds=(0, 1))

    alpha_init = {0: 1.95e-04, 1: 1e+06, 2: 2e+03, 3: 1e-03}

    #: Scalars
    m.cinit = Param(initialize=0.1367)
    m.tinit = Param(initialize=0.7293)
    m.uinit = Param(initialize=390.0)

    m.cdes = Param(initialize=0.0944)
    m.tdes = Param(initialize=0.7766)
    m.udes = Param(initialize=340)
    m.k10 = Param(initialize=300)
    m.n = Param(initialize=5)

    m.cf = Param(initialize=7.6)
    m.tf = Param(initialize=300)
    m.tc = Param(initialize=290)

    m.theta = Param(initialize=20)
    m.yf = Param(initialize=0.3947)
    m.yc = Param(initialize=0.3816)
    m.time = Param(initialize=10)
    m.point = Param(initialize=0)

    # m.nfe = Param(initialize=100)
    # m.ncp = Param(initialize=3)
    m.slopec = Param()
    m.slopet = Param()
    m.slopeu = Param()
    m.ii = Param()
    m.jj = Param()

    m.alpha = Param([0, 1, 2, 3], initialize=alpha_init)

    # m.cguess = Param(m.fe, m.cp)
    # m.tguess = Param(m.fe, m.cp)
    # m.ttguess = Param(m.fe, m.cp)
    # m.uguess = Param(m.fe, m.cp)

    point = 0
    slopec = (value(m.cdes) - value(m.cinit)) / (m.nfe * m.ncp)
    slopet = (value(m.tdes) - value(m.tinit)) / (m.nfe * m.ncp)
    slopeu = (value(m.udes) - value(m.uinit)) / (m.nfe * m.ncp)

    cguess = {}
    tguess = {}
    ttguess = {}
    uguess = {}

    for i in range(1, m.nfe + 1):
        for j in range(1, m.ncp + 1):
            point += 1
            cguess[i, j] = slopec * point + value(m.cinit)
            tguess[i, j] = slopet * point + value(m.tinit)
            ttguess[i, j] = value(m.time) * point
            uguess[i, j] = slopeu * point + value(m.uinit)

    #: One can pass the dictionary, as long as the keys are defined within the index set.
    m.C = Var(m.t, initialize=cguess)
    m.T = Var(m.t, initialize=tguess)
    m.u = Var(m.t, initialize=uguess)
    # m.tt = Var(m.t, initialize=ttguess)

    m.dC_dt = DerivativeVar(m.C)
    m.dT_dt = DerivativeVar(m.T)

    m.phi = Var()


    def _odec_rule(mod, t):
        if t > 0:
            return mod.dC_dt[t] == \
                   ((1 - mod.C[t]) / mod.theta - mod.k10 * exp(-mod.n / mod.T[t]) * mod.C[t])
        else:
            return Constraint.Skip


    def _odet_rule(mod, t):
        if t > 0:
            return mod.dT_dt[t] == \
                   ((mod.yf - mod.T[t]) / mod.theta + mod.k10 * exp(-mod.n / mod.T[t]) * mod.C[t] -
                    mod.alpha[0] * mod.u[t] * (mod.T[t] - mod.yc))
        else:
            return Constraint.Skip


    def _ic_rule(mod):
        return mod.C[0] == mod.cinit


    def _it_rule(mod):
        return mod.T[0] == mod.tinit


    m.OdeT = Constraint(m.t, rule=_odet_rule)
    m.OdeC = Constraint(m.t, rule=_odec_rule)

    m.IC = Constraint(rule=_ic_rule)
    m.IT = Constraint(rule=_it_rule)


    def objective_rule(mod):
        return sum((mod.alpha[1] * (mod.cdes - mod.C[t]) ** 2 +
                    mod.alpha[2] * (mod.tdes - mod.T[t]) ** 2 +
                    mod.alpha[3] * (mod.udes - mod.u[t]) ** 2)
                   for t in m.t if t > 0)


    m.fobj = Objective(sense=minimize, rule=objective_rule)

    dae = TransformationFactory('dae.collocation')
    dae.apply_to(m, nfe=m.nfe, ncp=m.ncp, scheme='LAGRANGE-RADAU')


    for var in m.C.itervalues():
        var.setlb(0)
        var.setub(1)

    for var in m.T.itervalues():
        var.setlb(0.1)
        var.setub(1)

    for var in m.u.itervalues():
        var.setlb(0)
        var.setub(500)

    for t in m.t:
        m.C[t].set_value(slopec * t + value(m.cinit))
        m.T[t].set_value(slopet * t + value(m.tinit))
        m.u[t].set_value(slopeu * t + value(m.uinit))

    ipopt = SolverFactory('ipopt')
    results = ipopt.solve(m, tee=True)
    #
    # if results.solver.status == SolverStatus.ok:
    #     print("Okay")
    #     ul = []
    #     tl = []
    #     for key in m.t:
    #         var = m.u[key]
    #         if var.stale:
    #             continue
    #         ul.append(value(var))
    #         tl.append(key)
    #     plt.plot(tl, ul)
    #     plt.show()

    if results.solver.status == SolverStatus.ok and with_plots:
        print("Okay")
        templ = []
        cl = []
        tl = []
        ul = []
        for key in m.t:
            var = m.u[key]
            if var.stale:
                continue
            templ.append(value(m.T[key]))
            cl.append(value(m.C[key]))
            tl.append(key)
            ul.append(value(m.u[key]))
        plt.subplot(3, 1, 1)
        plt.plot(tl, templ)
        plt.subplot(3, 1, 2)
        plt.plot(tl, cl)
        plt.subplot(3, 1, 3)
        plt.plot(tl, ul)
        plt.show()

if __name__ == '__main__':
    main()
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#

from __future__ import division
from __future__ import print_function

import matplotlib.pyplot as plt
from pyomo.core.kernel.expr import exp
from pyomo.core.kernel.numvalue import value
from pyomo.dae import *
#: pyomo imports
from pyomo.environ import *
from pyomo.opt import SolverFactory, SolverStatus

__author__ = "<NAME>"  #: May 2018

def main():
    with_plots = False
    #: Number of finite elements
    nfe = 100
    #: Number of collocation points
    ncp = 3
    m = ConcreteModel()

    m.nfe = nfe
    m.ncp = ncp

    m.t = ContinuousSet(bounds=(0, 1))

    alpha_init = {0: 1.95e-04, 1: 1e+06, 2: 2e+03, 3: 1e-03}

    #: Scalars
    m.cinit = Param(initialize=0.1367)
    m.tinit = Param(initialize=0.7293)
    m.uinit = Param(initialize=390.0)

    m.cdes = Param(initialize=0.0944)
    m.tdes = Param(initialize=0.7766)
    m.udes = Param(initialize=340)
    m.k10 = Param(initialize=300)
    m.n = Param(initialize=5)

    m.cf = Param(initialize=7.6)
    m.tf = Param(initialize=300)
    m.tc = Param(initialize=290)

    m.theta = Param(initialize=20)
    m.yf = Param(initialize=0.3947)
    m.yc = Param(initialize=0.3816)
    m.time = Param(initialize=10)
    m.point = Param(initialize=0)

    # m.nfe = Param(initialize=100)
    # m.ncp = Param(initialize=3)
    m.slopec = Param()
    m.slopet = Param()
    m.slopeu = Param()
    m.ii = Param()
    m.jj = Param()

    m.alpha = Param([0, 1, 2, 3], initialize=alpha_init)

    # m.cguess = Param(m.fe, m.cp)
    # m.tguess = Param(m.fe, m.cp)
    # m.ttguess = Param(m.fe, m.cp)
    # m.uguess = Param(m.fe, m.cp)

    point = 0
    slopec = (value(m.cdes) - value(m.cinit)) / (m.nfe * m.ncp)
    slopet = (value(m.tdes) - value(m.tinit)) / (m.nfe * m.ncp)
    slopeu = (value(m.udes) - value(m.uinit)) / (m.nfe * m.ncp)

    cguess = {}
    tguess = {}
    ttguess = {}
    uguess = {}

    for i in range(1, m.nfe + 1):
        for j in range(1, m.ncp + 1):
            point += 1
            cguess[i, j] = slopec * point + value(m.cinit)
            tguess[i, j] = slopet * point + value(m.tinit)
            ttguess[i, j] = value(m.time) * point
            uguess[i, j] = slopeu * point + value(m.uinit)

    #: One can pass the dictionary, as long as the keys are defined within the index set.
    m.C = Var(m.t, initialize=cguess)
    m.T = Var(m.t, initialize=tguess)
    m.u = Var(m.t, initialize=uguess)
    # m.tt = Var(m.t, initialize=ttguess)

    m.dC_dt = DerivativeVar(m.C)
    m.dT_dt = DerivativeVar(m.T)

    m.phi = Var()


    def _odec_rule(mod, t):
        if t > 0:
            return mod.dC_dt[t] == \
                   ((1 - mod.C[t]) / mod.theta - mod.k10 * exp(-mod.n / mod.T[t]) * mod.C[t])
        else:
            return Constraint.Skip


    def _odet_rule(mod, t):
        if t > 0:
            return mod.dT_dt[t] == \
                   ((mod.yf - mod.T[t]) / mod.theta + mod.k10 * exp(-mod.n / mod.T[t]) * mod.C[t] -
                    mod.alpha[0] * mod.u[t] * (mod.T[t] - mod.yc))
        else:
            return Constraint.Skip


    def _ic_rule(mod):
        return mod.C[0] == mod.cinit


    def _it_rule(mod):
        return mod.T[0] == mod.tinit


    m.OdeT = Constraint(m.t, rule=_odet_rule)
    m.OdeC = Constraint(m.t, rule=_odec_rule)

    m.IC = Constraint(rule=_ic_rule)
    m.IT = Constraint(rule=_it_rule)


    def objective_rule(mod):
        return sum((mod.alpha[1] * (mod.cdes - mod.C[t]) ** 2 +
                    mod.alpha[2] * (mod.tdes - mod.T[t]) ** 2 +
                    mod.alpha[3] * (mod.udes - mod.u[t]) ** 2)
                   for t in m.t if t > 0)


    m.fobj = Objective(sense=minimize, rule=objective_rule)

    dae = TransformationFactory('dae.collocation')
    dae.apply_to(m, nfe=m.nfe, ncp=m.ncp, scheme='LAGRANGE-RADAU')


    for var in m.C.itervalues():
        var.setlb(0)
        var.setub(1)

    for var in m.T.itervalues():
        var.setlb(0.1)
        var.setub(1)

    for var in m.u.itervalues():
        var.setlb(0)
        var.setub(500)

    for t in m.t:
        m.C[t].set_value(slopec * t + value(m.cinit))
        m.T[t].set_value(slopet * t + value(m.tinit))
        m.u[t].set_value(slopeu * t + value(m.uinit))

    ipopt = SolverFactory('ipopt')
    results = ipopt.solve(m, tee=True)
    #
    # if results.solver.status == SolverStatus.ok:
    #     print("Okay")
    #     ul = []
    #     tl = []
    #     for key in m.t:
    #         var = m.u[key]
    #         if var.stale:
    #             continue
    #         ul.append(value(var))
    #         tl.append(key)
    #     plt.plot(tl, ul)
    #     plt.show()

    if results.solver.status == SolverStatus.ok and with_plots:
        print("Okay")
        templ = []
        cl = []
        tl = []
        ul = []
        for key in m.t:
            var = m.u[key]
            if var.stale:
                continue
            templ.append(value(m.T[key]))
            cl.append(value(m.C[key]))
            tl.append(key)
            ul.append(value(m.u[key]))
        plt.subplot(3, 1, 1)
        plt.plot(tl, templ)
        plt.subplot(3, 1, 2)
        plt.plot(tl, cl)
        plt.subplot(3, 1, 3)
        plt.plot(tl, ul)
        plt.show()

if __name__ == '__main__':
    main()
#!/usr/bin/env python # -*- coding: utf-8 -*- # #: pyomo imports #: May 2018 #: Number of finite elements #: Number of collocation points #: Scalars # m.nfe = Param(initialize=100) # m.ncp = Param(initialize=3) # m.cguess = Param(m.fe, m.cp) # m.tguess = Param(m.fe, m.cp) # m.ttguess = Param(m.fe, m.cp) # m.uguess = Param(m.fe, m.cp) #: One can pass the dictionary, as long as the keys are defined within the index set. # m.tt = Var(m.t, initialize=ttguess) # # if results.solver.status == SolverStatus.ok: #     print("Okay") #     ul = [] #     tl = [] #     for key in m.t: #         var = m.u[key] #         if var.stale: #             continue #         ul.append(value(var)) #         tl.append(key) #     plt.plot(tl, ul) #     plt.show()