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

"""
***********************************************************************************
tutorial20.py
DAE Tools: pyDAE module, www.daetools.com
***********************************************************************************
DAE Tools is free software; you can redistribute it and/or modify it under the
Foundation. DAE Tools is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the
DAE Tools software; if not, see <http://www.gnu.org/licenses/>.
************************************************************************************
"""
__doc__ = """
This tutorial illustrates the support variable constraints available in Sundials IDA solver.
Benchmarks are available from `Matlab documentation
<https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html>`_.

1. Absolute Value Function:

dy/dt = -fabs(y)

solved on the interval [0,40] with the initial condition y(0) = 1.
The solution of this ODE decays to zero. If the solver produces a negative solution value,
the computation eventually will fail as the calculated solution diverges to -inf.
Using the constraint y >= 0 resolves this problem.

2. The Knee problem:

epsilon * dy/dt = (1-t)*y - y**2

solved on the interval [0,2] with the initial condition y(0) = 1.
The parameter epsilon is 0 < epsilon << 1 and in this example equal to 1e-6.
The solution follows the y = 1-x isocline for the whole interval of integration
which is incorrect. Using the constraint y >= 0 resolves the problem.

In DAE Tools contraints follow the Sundials IDA solver implementation and can be
specified using the valueConstraint argument of the daeVariableType class __init__
function:

- eNoConstraint (default)
- eValueGTEQ: imposes >= 0 constraint
- eValueLTEQ: imposes <= 0 constraint
- eValueGT:   imposes > 0 constraint
- eValueLT:   imposes < 0 constraint

and changed for individual variables using daeVariable.SetValueConstraint functions.

Absolute Value Function solution plot:

.. image:: _static/tutorial20-results1.png
:width: 500px

The Knee problem solution plot:

.. image:: _static/tutorial20-results2.png
:width: 500px
"""

import sys, numpy
from time import localtime, strftime
from daetools.pyDAE import *

# Impose >= constraint on y value using the eValueGTEQ flag.
type_y = daeVariableType("type_y", unit(), 0, 1E10, 0, 1e-5, eValueGTEQ)

class modTutorial1(daeModel):
def __init__(self, Name, Parent = None, Description = ""):
daeModel.__init__(self, Name, Parent, Description)

self.y = daeVariable("y", type_y, self, "")

def DeclareEquations(self):
daeModel.DeclareEquations(self)

# Auxiliary objects to make equations more readable
y     = self.y()
dy_dt = dt(self.y())

eq = self.CreateEquation("y")
eq.Residual = dy_dt + numpy.fabs(y)
eq.CheckUnitsConsistency = False

class modTutorial2(daeModel):
def __init__(self, Name, Parent = None, Description = ""):
daeModel.__init__(self, Name, Parent, Description)

self.y = daeVariable("y", type_y, self, "")

def DeclareEquations(self):
daeModel.DeclareEquations(self)

# 0 < eps << 1
epsilon = 1e-6

# Auxiliary objects to make equations more readable
t     = Time()
y     = self.y()
dy_dt = dt(self.y())

eq = self.CreateEquation("y")
eq.Residual = epsilon * dy_dt - ((1-t)*y - y**2)
eq.CheckUnitsConsistency = False

class simTutorial1(daeSimulation):
def __init__(self):
daeSimulation.__init__(self)
self.m = modTutorial1("tutorial20(1)")
self.m.Description = __doc__

def SetUpParametersAndDomains(self):
pass

def SetUpVariables(self):
self.m.y.SetInitialCondition(1.0)

class simTutorial2(daeSimulation):
def __init__(self):
daeSimulation.__init__(self)
self.m = modTutorial2("tutorial20(2)")
self.m.Description = __doc__

def SetUpParametersAndDomains(self):
pass

def SetUpVariables(self):
self.m.y.SetInitialCondition(1.0)

def run(**kwargs):
run1(**kwargs)
run2(**kwargs)

def run1(**kwargs):
simulation = simTutorial1()
return daeActivity.simulate(simulation, reportingInterval = 1.0,
timeHorizon       = 40.0,
**kwargs)

def run2(**kwargs):
simulation = simTutorial2()
return daeActivity.simulate(simulation, reportingInterval = 0.05,
timeHorizon       = 2.0,
**kwargs)

if __name__ == "__main__":
guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True
run(guiRun = guiRun)
```