00001 /* 00002 * This file is part of ACADO Toolkit. 00003 * 00004 * ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization. 00005 * Copyright (C) 2008-2009 by Boris Houska and Hans Joachim Ferreau, K.U.Leuven. 00006 * Developed within the Optimization in Engineering Center (OPTEC) under 00007 * supervision of Moritz Diehl. All rights reserved. 00008 * 00009 * ACADO Toolkit is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 3 of the License, or (at your option) any later version. 00013 * 00014 * ACADO Toolkit is distributed in the hope that it will be useful, 00015 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 * Lesser General Public License for more details. 00018 * 00019 * You should have received a copy of the GNU Lesser General Public 00020 * License along with ACADO Toolkit; if not, write to the Free Software 00021 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00022 * 00023 */ 00024 00025 00034 void glycemia( DifferentialEquation *f ){ 00035 00036 // Define a Right-Hand-Side: 00037 // ------------------------- 00038 00039 DifferentialState G ; // blood glucose concentration 00040 DifferentialState I ; // blood insulin concentration 00041 DifferentialState X ; // effect of insulin on net glucose disappearance 00042 DifferentialState I2; // effect of endogenous insulin 00043 00044 Control FI; // exogenous insulin flow 00045 00046 Disturbance w0; 00047 Disturbance w1; 00048 Disturbance w2; 00049 Disturbance w3; 00050 00051 Parameter FG; // carbohydrate calories flow 00052 Parameter P1; 00053 Parameter P2; 00054 Parameter P3; 00055 Parameter C1; 00056 Parameter C2; 00057 Parameter C3; 00058 Parameter C4; 00059 Parameter C5; 00060 Parameter n; 00061 Parameter alpha; 00062 Parameter gamma; 00063 00064 const double s = 0.1 ; 00065 00066 IntermediateState z; 00067 00068 z = exp(I2/s); 00069 00070 *f << dot(G) == (-P1 - X) * G + P1 * C1 + FG/C2 + w0; 00071 *f << dot(I) == alpha * s *log(1 + z) - n*(I-C4) + FI/C3 + w1; 00072 *f << dot(X) == -P2 * X + P3 * 0.001 * (I-C4) + w2; 00073 *f << dot(I2) == gamma * (G - C5) - n * (I2) + w3; 00074 } 00075