dirk_export.cpp
Go to the documentation of this file.
1 /*
2  * This file is part of ACADO Toolkit.
3  *
4  * ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
5  * Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
6  * Milan Vukov, Rien Quirynen, KU Leuven.
7  * Developed within the Optimization in Engineering Center (OPTEC)
8  * under supervision of Moritz Diehl. All rights reserved.
9  *
10  * ACADO Toolkit is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 3 of the License, or (at your option) any later version.
14  *
15  * ACADO Toolkit is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  * Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with ACADO Toolkit; if not, write to the Free Software
22  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23  *
24  */
25 
26 
27 
35 
36 #include <sstream>
37 using namespace std;
38 
39 
40 
42 
43 
44 //
45 // PUBLIC MEMBER FUNCTIONS:
46 //
47 
49  const std::string& _commonHeaderName
50  ) : ForwardIRKExport( _userInteraction,_commonHeaderName )
51 {
52 
53 }
54 
56 {
57 
58 }
59 
60 
62 {
63  if ( solver )
64  delete solver;
65  solver = 0;
66 
67  clear( );
68 }
69 
70 
72 
73  if( this != &arg ){
74 
76  }
77  return *this;
78 }
79 
80 
82 {
83  if( NX1 > 0 ) {
84  ExportForLoop loop( index1,0,numStages );
85  loop.addStatement( rk_xxx.getCols(0,NX1) == rk_eta.getCols(0,NX1) );
86  ExportForLoop loop01( index2,0,NX1 );
87  ExportForLoop loop02( index3,0,index1 );
88  loop02.addStatement( rk_xxx.getCol( index2 ) += Ah.getElement(index1,index3)*rk_kkk.getElement(index2,index3) );
89  loop01.addStatement( loop02 );
90  loop.addStatement( loop01 );
92 
93  ExportForLoop loop5( index2,0,NX1 );
94  loop5.addStatement( tmp_index == index1*NX1+index2 );
95  loop5.addStatement( rk_kkk.getElement(index2,index1) == rk_mat1.getElement(tmp_index,0)*rk_b.getRow(0) );
96  ExportForLoop loop6( index3,1,NX1 );
97  loop6.addStatement( rk_kkk.getElement(index2,index1) += rk_mat1.getElement(tmp_index,index3)*rk_b.getRow(index3) );
98  loop5.addStatement(loop6);
99  loop.addStatement(loop5);
100  block->addStatement(loop);
101  }
102 
103  return SUCCESSFUL_RETURN;
104 }
105 
106 
108 {
109  if( NX1 > 0 ) {
110  DMatrix mat1 = formMatrix( M11, A11 );
111  rk_mat1 = ExportVariable( "rk_mat1", mat1, STATIC_CONST_REAL );
112  code.addDeclaration( rk_mat1 );
113  // TODO: Ask Milan why this does NOT work properly !!
115  double h = (grid.getLastTime() - grid.getFirstTime())/grid.getNumIntervals();
116 
117  DMatrix sens = zeros<double>(NX1*(NX1+NU), numStages);
118  uint i, j, k, s1, s2;
119  for( i = 0; i < NX1; i++ ) {
120  DVector vec(NX1);
121  for( j = 0; j < numStages; j++ ) {
122  for( k = 0; k < NX1; k++ ) {
123  vec(k) = A11(k,i);
124  for( s1 = 0; s1 < j; s1++ ) {
125  for( s2 = 0; s2 < NX1; s2++ ) {
126  vec(k) = vec(k) + AA(j,s1)*h*A11(k,s2)*sens(i*NX1+s2,s1);
127  }
128  }
129  }
130  DVector sol = mat1*vec;
131  for( k = 0; k < NX1; k++ ) {
132  sens(i*NX1+k,j) = sol(k);
133  }
134  }
135  }
136  for( i = 0; i < NU; i++ ) {
137  DVector vec(NX1);
138  for( j = 0; j < numStages; j++ ) {
139  for( k = 0; k < NX1; k++ ) {
140  vec(k) = B11(k,i);
141  for( s1 = 0; s1 < j; s1++ ) {
142  for( s2 = 0; s2 < NX1; s2++ ) {
143  vec(k) = vec(k) + AA(j,s1)*h*A11(k,s2)*sens(NX1*NX1+i*NX1+s2,s1);
144  }
145  }
146  }
147  DVector sol = mat1*vec;
148  for( k = 0; k < NX1; k++ ) {
149  sens(NX1*NX1+i*NX1+k,j) = sol(k);
150  }
151  }
152  }
153  rk_dk1 = ExportVariable( "rk_dk1", sens, STATIC_CONST_REAL );
154  code.addDeclaration( rk_dk1 );
155  // TODO: Ask Milan why this does NOT work properly !!
156  rk_dk1 = ExportVariable( "rk_dk1", NX1*(NX1+NU), numStages, STATIC_CONST_REAL, ACADO_LOCAL );
157  }
158 
159  return SUCCESSFUL_RETURN;
160 }
161 
162 
164  if( jacobian.getNumRows() != jacobian.getNumCols() ) {
166  }
167  double h = (grid.getLastTime() - grid.getFirstTime())/grid.getNumIntervals();
168  uint vars = jacobian.getNumRows();
169  uint i1, i2, j2;
170  DMatrix result = zeros<double>(numStages*vars, vars);
171  DMatrix tmp = zeros<double>(vars, vars);
172  for( i1 = 0; i1 < numStages; i1++ ){
173  for( i2 = 0; i2 < vars; i2++ ){
174  for( j2 = 0; j2 < vars; j2++ ) {
175  tmp(i2, j2) = mass(i2,j2) - AA(i1,i1)*h*jacobian(i2,j2);
176  }
177  }
178  tmp = tmp.inverse();
179  for( i2 = 0; i2 < vars; i2++ ){
180  for( j2 = 0; j2 < vars; j2++ ) {
181  result(i1*vars+i2, j2) = tmp(i2, j2);
182  }
183  }
184  }
185 
186  return result;
187 }
188 
189 
190 returnValue DiagonallyImplicitRKExport::solveImplicitSystem( ExportStatementBlock* block, const ExportIndex& index1, const ExportIndex& index2, const ExportIndex& index3, const ExportIndex& tmp_index, const ExportVariable& Ah, const ExportVariable& C, const ExportVariable& det, bool DERIVATIVES )
191 {
192  if( NX2 > 0 || NXA > 0 ) {
193 
194  if( REUSE ) block->addStatement( std::string( "if( " ) + reset_int.getFullName() + " ) {\n" );
195  // Initialization iterations:
196  ExportForLoop loop11( index2,0,numStages );
197  ExportForLoop loop1( index1,0,numItsInit+1 ); // NOTE: +1 because 0 will lead to NaNs, so the minimum number of iterations is 1 at the initialization
198  evaluateMatrix( &loop1, index2, index3, tmp_index, rk_A, Ah, C, true, DERIVATIVES );
199  loop1.addStatement( det.getFullName() + " = " + ExportStatement::fcnPrefix + "_" + solver->getNameSolveFunction() + "( &" + rk_A.get(index2*(NX2+NXA),0) + ", " + rk_b.getFullName() + ", &" + rk_auxSolver.get(index2,0) + " );\n" );
200  loop1.addStatement( rk_kkk.getSubMatrix( NX1,NX1+NX2,index2,index2+1 ) += rk_b.getRows( 0,NX2 ) ); // differential states
201  if(NXA > 0) loop1.addStatement( rk_kkk.getSubMatrix( NX,NX+NXA,index2,index2+1 ) += rk_b.getRows( NX2,NX2+NXA ) ); // algebraic states
202  loop11.addStatement( loop1 );
203  block->addStatement( loop11 );
204  if( REUSE ) block->addStatement( std::string( "}\n" ) );
205 
206  // the rest (numIts) of the Newton iterations with reuse of the Jacobian (no evaluation or factorization needed)
207  ExportForLoop loop21( index2,0,numStages );
208  ExportForLoop loop2( index1,0,numIts );
209  evaluateStatesImplicitSystem( &loop2, Ah, C, index2, index3, tmp_index );
210  evaluateRhsImplicitSystem( &loop2, index2 );
212  loop2.addStatement( rk_kkk.getSubMatrix( NX1,NX1+NX2,index2,index2+1 ) += rk_b.getRows( 0,NX2 ) ); // differential states
213  if(NXA > 0) loop2.addStatement( rk_kkk.getSubMatrix( NX,NX+NXA,index2,index2+1 ) += rk_b.getRows( NX2,NX2+NXA ) ); // algebraic states
214  loop21.addStatement( loop2 );
215  block->addStatement( loop21 );
216 
217  if( DERIVATIVES ) {
218  // solution calculated --> evaluate and save the necessary derivatives in rk_diffsTemp and update the matrix rk_A:
219  ExportForLoop loop3( index2,0,numStages );
220  evaluateMatrix( &loop3, index2, index3, tmp_index, rk_A, Ah, C, false, DERIVATIVES );
221  block->addStatement( loop3 );
222  }
223 
224  // IF DEBUG MODE:
225  int debugMode;
226  get( INTEGRATOR_DEBUG_MODE, debugMode );
227  if ( (bool)debugMode == true ) {
228  block->addStatement( debug_mat == rk_A );
229  }
230  }
231 
232  return SUCCESSFUL_RETURN;
233 }
234 
235 
236 returnValue DiagonallyImplicitRKExport::sensitivitiesImplicitSystem( ExportStatementBlock* block, const ExportIndex& index1, const ExportIndex& index2, const ExportIndex& index3, const ExportIndex& tmp_index1, const ExportIndex& tmp_index2, const ExportVariable& Ah, const ExportVariable& Bh, const ExportVariable& det, bool STATES, uint number )
237 {
238  if( NX2 > 0 ) {
239  DMatrix zeroM = zeros<double>( NX2+NXA,1 );
240  DMatrix tempCoefs( evaluateDerivedPolynomial( 0.0 ) );
241  uint i;
242 
243  ExportForLoop loop1( index2,0,numStages );
244  if( STATES && number == 1 ) {
245  ExportForLoop loop2( index3,0,NX1 );
246  loop2.addStatement( std::string(rk_rhsTemp.get( index3,0 )) + " = -(" + index3.getName() + " == " + index1.getName() + ");\n" );
247  ExportForLoop loop21( tmp_index1,0,index2+1 );
248  loop21.addStatement( rk_rhsTemp.getRow( index3 ) -= rk_diffK.getElement( index3,tmp_index1 )*Ah.getElement(index2,tmp_index1) );
249  loop2.addStatement( loop21 );
250  loop1.addStatement( loop2 );
251  ExportForLoop loop3( index3,0,NX2+NXA );
252  loop3.addStatement( rk_b.getRow( index3 ) == rk_diffsTemp2.getSubMatrix( index2,index2+1,index3*(NVARS2),index3*(NVARS2)+NX1 )*rk_rhsTemp.getRows(0,NX1) );
253  if( NDX2 > 0 ) {
254  loop3.addStatement( rk_b.getRow( index3 ) -= rk_diffsTemp2.getSubMatrix( index2,index2+1,index3*(NVARS2)+NVARS2-NX1-NX2,index3*(NVARS2)+NVARS2-NX2 )*rk_diffK.getSubMatrix( 0,NX1,index2,index2+1 ) );
255  }
256  loop1.addStatement( loop3 );
257  }
258  else if( STATES && number == 2 ) {
259  for( i = 0; i < NX2+NXA; i++ ) {
260  loop1.addStatement( rk_b.getRow( i ) == zeroM.getRow( 0 ) - rk_diffsTemp2.getElement( index2,index1+i*(NVARS2) ) );
261  }
262  }
263  else { // ~= STATES
264  ExportForLoop loop2( index3,0,NX1 );
265  loop2.addStatement( rk_rhsTemp.getRow( index3 ) == rk_diffK.getElement( index3,0 )*Ah.getElement(index2,0) );
266  ExportForLoop loop21( tmp_index1,1,index2+1 );
267  loop21.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,tmp_index1 )*Ah.getElement(index2,tmp_index1) );
268  loop2.addStatement( loop21 );
269  loop1.addStatement( loop2 );
270  ExportForLoop loop3( index3,0,NX2+NXA );
271  loop3.addStatement( tmp_index2 == index1+index3*(NVARS2) );
272  loop3.addStatement( rk_b.getRow( index3 ) == zeroM.getRow( 0 ) - rk_diffsTemp2.getElement( index2,tmp_index2+NX1+NX2+NXA ) );
273  loop3.addStatement( rk_b.getRow( index3 ) -= rk_diffsTemp2.getSubMatrix( index2,index2+1,index3*(NVARS2),index3*(NVARS2)+NX1 )*rk_rhsTemp.getRows(0,NX1) );
274  if( NDX2 > 0 ) {
275  loop3.addStatement( rk_b.getRow( index3 ) -= rk_diffsTemp2.getSubMatrix( index2,index2+1,index3*(NVARS2)+NVARS2-NX1-NX2,index3*(NVARS2)+NVARS2-NX2 )*rk_diffK.getSubMatrix( 0,NX1,index2,index2+1 ) );
276  }
277  loop1.addStatement( loop3 );
278  }
279  ExportForLoop loop11( index3,0,NX2+NXA );
280  ExportForLoop loop12( tmp_index1,0,index2 );
281  ExportForLoop loop13( tmp_index2,NX1,NX1+NX2 );
282  loop13.addStatement( std::string( rk_b.get(index3,0) ) + " -= " + Ah.get(index2,tmp_index1) + "*" + rk_diffsTemp2.get(index2,index3*NVARS2+tmp_index2) + "*" + rk_diffK.get(tmp_index2,tmp_index1) + ";\n" );
283  loop12.addStatement( loop13 );
284  loop11.addStatement( loop12 );
285  loop1.addStatement( loop11 );
286  if( STATES && (number == 1 || NX1 == 0) ) {
287  loop1.addStatement( std::string( "if( 0 == " ) + index1.getName() + " ) {\n" ); // factorization of the new matrix rk_A not yet calculated!
288  loop1.addStatement( det.getFullName() + " = " + ExportStatement::fcnPrefix + "_" + solver->getNameSolveFunction() + "( &" + rk_A.get(index2*(NX2+NXA),0) + ", " + rk_b.getFullName() + ", &" + rk_auxSolver.get(index2,0) + " );\n" );
289  loop1.addStatement( std::string( "}\n else {\n" ) );
290  }
292  if( STATES && (number == 1 || NX1 == 0) ) loop1.addStatement( std::string( "}\n" ) );
293  // update rk_diffK with the new sensitivities:
294  loop1.addStatement( rk_diffK.getSubMatrix(NX1,NX1+NX2,index2,index2+1) == rk_b.getRows(0,NX2) );
295  loop1.addStatement( rk_diffK.getSubMatrix(NX,NX+NXA,index2,index2+1) == rk_b.getRows(NX2,NX2+NXA) );
296  block->addStatement( loop1 );
297  // update rk_diffsNew with the new sensitivities:
298  ExportForLoop loop3( index2,0,NX2 );
299  if( STATES && number == 2 ) loop3.addStatement( std::string(rk_diffsNew2.get( index2,index1 )) + " = (" + index2.getName() + " == " + index1.getName() + "-" + toString(NX1) + ");\n" );
300 
301  if( STATES && number == 2 ) loop3.addStatement( rk_diffsNew2.getElement( index2,index1 ) += rk_diffK.getRow( NX1+index2 )*Bh );
302  else if( STATES ) loop3.addStatement( rk_diffsNew2.getElement( index2,index1 ) == rk_diffK.getRow( NX1+index2 )*Bh );
303  else loop3.addStatement( rk_diffsNew2.getElement( index2,index1+NX1+NX2 ) == rk_diffK.getRow( NX1+index2 )*Bh );
304  block->addStatement( loop3 );
305  if( NXA > 0 ) {
306  block->addStatement( std::string("if( run == 0 ) {\n") );
307  ExportForLoop loop4( index2,0,NXA );
308  if( STATES ) loop4.addStatement( rk_diffsNew2.getElement( index2+NX2,index1 ) == rk_diffK.getRow( NX+index2 )*tempCoefs );
309  else loop4.addStatement( rk_diffsNew2.getElement( index2+NX2,index1+NX1+NX2 ) == rk_diffK.getRow( NX+index2 )*tempCoefs );
310  block->addStatement( loop4 );
311  block->addStatement( std::string("}\n") );
312  }
313  }
314 
315  return SUCCESSFUL_RETURN;
316 }
317 
318 
319 returnValue DiagonallyImplicitRKExport::evaluateMatrix( ExportStatementBlock* block, const ExportIndex& index1, const ExportIndex& index2, const ExportIndex& tmp_index, const ExportVariable& _rk_A, const ExportVariable& Ah, const ExportVariable& C, bool evaluateB, bool DERIVATIVES )
320 {
321  evaluateStatesImplicitSystem( block, Ah, C, index1, index2, tmp_index );
322 
323  ExportIndex indexDiffs(index1);
324  if( !DERIVATIVES ) indexDiffs = ExportIndex(0);
325 
326  block->addFunctionCall( getNameDiffsRHS(), rk_xxx, rk_diffsTemp2.getAddress(indexDiffs,0) );
327  ExportForLoop loop2( index2,0,NX2+NXA );
328  loop2.addStatement( tmp_index == index1*(NX2+NXA)+index2 );
329  if( NDX2 == 0 ) {
330  loop2.addStatement( _rk_A.getSubMatrix( tmp_index,tmp_index+1,0,NX2 ) == Ah.getElement( 0,0 )*rk_diffsTemp2.getSubMatrix( indexDiffs,indexDiffs+1,index2*(NVARS2)+NX1,index2*(NVARS2)+NX1+NX2 ) );
331  loop2.addStatement( _rk_A.getElement( tmp_index,index2 ) -= 1 );
332  }
333  else {
334  loop2.addStatement( _rk_A.getSubMatrix( tmp_index,tmp_index+1,0,NX2 ) == Ah.getElement( 0,0 )*rk_diffsTemp2.getSubMatrix( indexDiffs,indexDiffs+1,index2*(NVARS2)+NX1,index2*(NVARS2)+NX1+NX2 ) );
335  loop2.addStatement( _rk_A.getSubMatrix( tmp_index,tmp_index+1,0,NX2 ) += rk_diffsTemp2.getSubMatrix( indexDiffs,indexDiffs+1,index2*(NVARS2)+NVARS2-NX2,index2*(NVARS2)+NVARS2 ) );
336  }
337  if( NXA > 0 ) {
338  DMatrix zeroM = zeros<double>( 1,NXA );
339  loop2.addStatement( _rk_A.getSubMatrix( tmp_index,tmp_index+1,NX2,NX2+NXA ) == rk_diffsTemp2.getSubMatrix( indexDiffs,indexDiffs+1,index2*(NVARS2)+NX1+NX2,index2*(NVARS2)+NX1+NX2+NXA ) );
340  }
341  block->addStatement( loop2 );
342  if( evaluateB ) {
343  evaluateRhsImplicitSystem( block, index1 );
344  }
345 
346  return SUCCESSFUL_RETURN;
347 }
348 
349 
351 {
352  ExportForLoop loop1( i, 0, NX1+NX2 );
353  loop1.addStatement( rk_xxx.getCol( i ) == rk_eta.getCol( i ) );
354  ExportForLoop loop2( j, 0, stage+1 );
355  loop2.addStatement( rk_xxx.getCol( i ) += Ah.getElement(stage,j)*rk_kkk.getElement( i,j ) );
356  loop1.addStatement( loop2 );
357  block->addStatement( loop1 );
358 
359  ExportForLoop loop3( i, 0, NXA );
360  loop3.addStatement( rk_xxx.getCol( NX+i ) == rk_kkk.getElement( NX+i,stage ) );
361  block->addStatement( loop3 );
362 
363  ExportForLoop loop4( i, 0, NDX2 );
364  loop4.addStatement( rk_xxx.getCol( inputDim-diffsDim+i ) == rk_kkk.getElement( i,stage ) );
365  block->addStatement( loop4 );
366 
367  if( C.getDim() > 0 ) { // There is a time dependence, so it must be set
368  block->addStatement( rk_xxx.getCol( inputDim-diffsDim+NDX2 ) == C.getCol(stage) );
369  }
370 
371  return SUCCESSFUL_RETURN;
372 }
373 
374 
376 {
377  DMatrix zeroM = zeros<double>( NX2+NXA,1 );
379  // matrix rk_b:
380  if( NDX2 == 0 ) {
381  block->addStatement( rk_b.getRows( 0,NX2 ) == rk_kkk.getSubMatrix( NX1,NX1+NX2,stage,stage+1 ) - rk_rhsTemp.getRows( 0,NX2 ) );
382  }
383  else {
384  block->addStatement( rk_b.getRows( 0,NX2 ) == zeroM.getRows( 0,NX2-1 ) - rk_rhsTemp.getRows( 0,NX2 ) );
385  }
386  if( NXA > 0 ) {
387  block->addStatement( rk_b.getRows( NX2,NX2+NXA ) == zeroM.getRows( 0,NXA-1 ) - rk_rhsTemp.getRows( NX2,NX2+NXA ) );
388  }
389 
390  return SUCCESSFUL_RETURN;
391 }
392 
393 
394 returnValue DiagonallyImplicitRKExport::solveOutputSystem( ExportStatementBlock* block, const ExportIndex& index1, const ExportIndex& index2, const ExportIndex& index3, const ExportIndex& tmp_index, const ExportVariable& Ah, bool DERIVATIVES )
395 {
396  if( NX3 > 0 ) {
397  ExportForLoop loop( index1,0,numStages );
398  evaluateStatesOutputSystem( &loop, Ah, index1 );
399  ExportForLoop loop01( index2,NX1+NX2,NX );
400  ExportForLoop loop02( index3,0,index1 );
401  loop02.addStatement( rk_xxx.getCol( index2 ) += Ah.getElement(index1,index3)*rk_kkk.getElement(index2,index3) );
402  loop01.addStatement( loop02 );
403  loop.addStatement( loop01 );
405  if( DERIVATIVES ) loop.addFunctionCall( getNameOutputDiffs(), rk_xxx, rk_diffsTemp3.getAddress(index1,0) );
406 
407  ExportForLoop loop5( index2,0,NX3 );
408  loop5.addStatement( tmp_index == index1*NX3+index2 );
409  loop5.addStatement( rk_kkk.getElement(NX1+NX2+index2,index1) == rk_mat3.getElement(tmp_index,0)*rk_b.getRow(0) );
410  ExportForLoop loop6( index3,1,NX3 );
411  loop6.addStatement( rk_kkk.getElement(NX1+NX2+index2,index1) += rk_mat3.getElement(tmp_index,index3)*rk_b.getRow(index3) );
412  loop5.addStatement(loop6);
413  loop.addStatement(loop5);
414  block->addStatement(loop);
415  }
416 
417  return SUCCESSFUL_RETURN;
418 }
419 
420 
421 returnValue DiagonallyImplicitRKExport::sensitivitiesOutputSystem( ExportStatementBlock* block, const ExportIndex& index1, const ExportIndex& index2, const ExportIndex& index3, const ExportIndex& index4, const ExportIndex& tmp_index1, const ExportIndex& tmp_index2, const ExportVariable& Ah, const ExportVariable& Bh, bool STATES, uint number )
422 {
423  if( NX3 > 0 ) {
424  uint i, j;
425  ExportForLoop loop1( index2,0,numStages );
426  if( STATES && number == 1 ) {
427  ExportForLoop loop2( index3,0,NX1 );
428  loop2.addStatement( std::string(rk_rhsTemp.get( index3,0 )) + " = (" + index3.getName() + " == " + index1.getName() + ");\n" );
429  for( i = 0; i < numStages; i++ ) {
430  loop2.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,i )*Ah.getElement(index2,i) );
431  }
432  loop1.addStatement( loop2 );
433  ExportForLoop loop3( index3,NX1,NX1+NX2 );
434  loop3.addStatement( rk_rhsTemp.getRow( index3 ) == 0.0 );
435  for( i = 0; i < numStages; i++ ) {
436  loop3.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,i )*Ah.getElement(index2,i) );
437  }
438  loop1.addStatement( loop3 );
439  ExportForLoop loop4( index3,0,NX3 );
440  loop4.addStatement( rk_b.getRow( index3 ) == rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3),index3*(NVARS3)+NX1+NX2 )*rk_rhsTemp.getRows(0,NX1+NX2) );
441  if( NXA3 > 0 ) {
442  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NX1+NX2,index3*(NVARS3)+NX1+NX2+NXA )*rk_diffK.getSubMatrix( NX,NX+NXA,index2,index2+1 ) );
443  }
444  if( NDX3 > 0 ) {
445  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NVARS3-NX1-NX2,index3*(NVARS3)+NVARS3 )*rk_diffK.getSubMatrix( 0,NX1+NX2,index2,index2+1 ) );
446  }
447  loop1.addStatement( loop4 );
448  }
449  else if( STATES && number == 2 ) {
450  ExportForLoop loop3( index3,NX1,NX1+NX2 );
451  loop3.addStatement( std::string(rk_rhsTemp.get( index3,0 )) + " = (" + index3.getName() + " == " + index1.getName() + ");\n" );
452  for( i = 0; i < numStages; i++ ) {
453  loop3.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,i )*Ah.getElement(index2,i) );
454  }
455  loop1.addStatement( loop3 );
456  ExportForLoop loop4( index3,0,NX3 );
457  loop4.addStatement( rk_b.getRow( index3 ) == rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NX1,index3*(NVARS3)+NX1+NX2 )*rk_rhsTemp.getRows(NX1,NX1+NX2) );
458  if( NXA3 > 0 ) {
459  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NX1+NX2,index3*(NVARS3)+NX1+NX2+NXA )*rk_diffK.getSubMatrix( NX,NX+NXA,index2,index2+1 ) );
460  }
461  if( NDX3 > 0 ) {
462  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NVARS3-NX2,index3*(NVARS3)+NVARS3 )*rk_diffK.getSubMatrix( NX1,NX1+NX2,index2,index2+1 ) );
463  }
464  loop1.addStatement( loop4 );
465  }
466  else if( !STATES ) {
467  ExportForLoop loop2( index3,0,NX1 );
468  loop2.addStatement( rk_rhsTemp.getRow( index3 ) == rk_diffK.getElement( index3,0 )*Ah.getElement(index2,0) );
469  for( i = 1; i < numStages; i++ ) {
470  loop2.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,i )*Ah.getElement(index2,i) );
471  }
472  loop1.addStatement( loop2 );
473  ExportForLoop loop3( index3,NX1,NX1+NX2 );
474  loop3.addStatement( rk_rhsTemp.getRow( index3 ) == rk_diffK.getElement( index3,0 )*Ah.getElement(index2,0) );
475  for( i = 1; i < numStages; i++ ) {
476  loop3.addStatement( rk_rhsTemp.getRow( index3 ) += rk_diffK.getElement( index3,i )*Ah.getElement(index2,i) );
477  }
478  loop1.addStatement( loop3 );
479  ExportForLoop loop4( index3,0,NX3 );
480  loop4.addStatement( tmp_index2 == index1+index3*(NVARS3) );
481  loop4.addStatement( rk_b.getRow( index3 ) == rk_diffsTemp3.getElement( index2,tmp_index2+NX1+NX2+NXA3 ) );
482  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3),index3*(NVARS3)+NX1+NX2 )*rk_rhsTemp.getRows(0,NX1+NX2) );
483  if( NXA3 > 0 ) {
484  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NX1+NX2,index3*(NVARS3)+NX1+NX2+NXA )*rk_diffK.getSubMatrix( NX,NX+NXA,index2,index2+1 ) );
485  }
486  if( NDX3 > 0 ) {
487  loop4.addStatement( rk_b.getRow( index3 ) += rk_diffsTemp3.getSubMatrix( index2,index2+1,index3*(NVARS3)+NVARS3-NX1-NX2,index3*(NVARS3)+NVARS3 )*rk_diffK.getSubMatrix( 0,NX1+NX2,index2,index2+1 ) );
488  }
489  loop1.addStatement( loop4 );
490  }
491  if( !STATES || number != 3 ) {
492  ExportForLoop loop12( tmp_index1,0,index2 );
493  for( i = 0; i < NX3; i++ ) {
494  for( j = NX1+NX2; j < NX; j++ ) {
495  if( acadoRoundAway(A33(i,j-NX1-NX2)) != 0 ) {
496  loop12.addStatement( std::string( rk_b.get(i,0) ) + " += " + Ah.get(index2,tmp_index1) + "*" + toString(A33(i,j-NX1-NX2)) + "*" + rk_diffK.get(j,tmp_index1) + ";\n" );
497  }
498  }
499  }
500  loop1.addStatement( loop12 );
501  }
502 
503  // update rk_diffK with the new sensitivities:
504  if( STATES && number == 3 ) {
505  block->addStatement( rk_diffK.getRows(NX1+NX2,NX) == rk_dk3.getRows(index1*NX3-(NX1+NX2)*NX3,index1*NX3+NX3-(NX1+NX2)*NX3) );
506  }
507  else {
508  ExportForLoop loop5( index3,0,NX3 );
509  loop5.addStatement( tmp_index1 == index2*NX3+index3 );
510  loop5.addStatement( rk_diffK.getElement(NX1+NX2+index3,index2) == rk_mat3.getElement(tmp_index1,0)*rk_b.getRow(0) );
511  ExportForLoop loop6( index4,1,NX3 );
512  loop6.addStatement( rk_diffK.getElement(NX1+NX2+index3,index2) += rk_mat3.getElement(tmp_index1,index4)*rk_b.getRow(index4) );
513  loop5.addStatement(loop6);
514  loop1.addStatement(loop5);
515  block->addStatement( loop1 );
516  }
517  // update rk_diffsNew with the new sensitivities:
518  if( grid.getNumIntervals() > 1 || !equidistantControlGrid() ) block->addStatement( std::string( "if( run == 0 ) {\n" ) );
519  ExportForLoop loop8( index2,0,NX3 );
520  if( STATES && number == 3 ) loop8.addStatement( std::string(rk_diffsNew3.get( index2,index1 )) + " = (" + index2.getName() + " == " + index1.getName() + "-" + toString(NX1+NX2) + ");\n" );
521 
522  if( STATES && number == 3 ) loop8.addStatement( rk_diffsNew3.getElement( index2,index1 ) += rk_diffK.getRow( NX1+NX2+index2 )*Bh );
523  else if( STATES ) loop8.addStatement( rk_diffsNew3.getElement( index2,index1 ) == rk_diffK.getRow( NX1+NX2+index2 )*Bh );
524  else loop8.addStatement( rk_diffsNew3.getElement( index2,index1+NX ) == rk_diffK.getRow( NX1+NX2+index2 )*Bh );
525  block->addStatement( loop8 );
526  if( grid.getNumIntervals() > 1 || !equidistantControlGrid() ) block->addStatement( std::string( "}\n" ) );
527  }
528 
529  return SUCCESSFUL_RETURN;
530 }
531 
532 
534 {
535  if( NX3 > 0 ) {
536  DMatrix mat3 = formMatrix( M33, A33 );
537  rk_mat3 = ExportVariable( "rk_mat3", mat3, STATIC_CONST_REAL );
538  code.addDeclaration( rk_mat3 );
539  // TODO: Ask Milan why this does NOT work properly !!
541  double h = (grid.getLastTime() - grid.getFirstTime())/grid.getNumIntervals();
542 
543  DMatrix sens = zeros<double>(NX3*NX3, numStages);
544  uint i, j, k, s1, s2;
545  for( i = 0; i < NX3; i++ ) {
546  DVector vec(NX3);
547  for( j = 0; j < numStages; j++ ) {
548  for( k = 0; k < NX3; k++ ) {
549  vec(k) = A33(k,i);
550  for( s1 = 0; s1 < j; s1++ ) {
551  for( s2 = 0; s2 < NX3; s2++ ) {
552  vec(k) = vec(k) + AA(j,s1)*h*A33(k,s2)*sens(i*NX3+s2,s1);
553  }
554  }
555  }
556  DVector sol = mat3*vec;
557  for( k = 0; k < NX3; k++ ) {
558  sens(i*NX3+k,j) = sol(k);
559  }
560  }
561  }
562  rk_dk3 = ExportVariable( "rk_dk3", sens, STATIC_CONST_REAL );
563  code.addDeclaration( rk_dk3 );
564  // TODO: Ask Milan why this does NOT work properly !!
565  rk_dk3 = ExportVariable( "rk_dk3", NX3*NX3, numStages, STATIC_CONST_REAL, ACADO_LOCAL );
566  }
567 
568  return SUCCESSFUL_RETURN;
569 }
570 
571 
573 {
575 
576  int debugMode;
577  get( INTEGRATOR_DEBUG_MODE, debugMode );
578 
579  int useOMP;
580  get(CG_USE_OPENMP, useOMP);
581  ExportStruct structWspace;
582  structWspace = useOMP ? ACADO_LOCAL : ACADO_WORKSPACE;
583 
584  rk_A = ExportVariable( "rk_A", numStages*(NX2+NXA), NX2+NXA, REAL, structWspace );
585  if ( (bool)debugMode == true && useOMP ) {
586  return ACADOERROR( RET_INVALID_OPTION );
587  }
588  else {
590  }
591  uint Xmax = NX1;
592  if( NX2 > Xmax ) Xmax = NX2;
593  if( NX3 > Xmax ) Xmax = NX3;
594  rk_b = ExportVariable( "rk_b", Xmax+NXA, 1, REAL, structWspace );
595 
596  if( NX2 > 0 || NXA > 0 ) {
597  solver->init( NX2+NXA );
598  solver->setup();
600  }
601 
602  return IRKsetup;
603 }
604 
605 
607 
608 // end of file.
ExportVariable debug_mat
Definition: irk_export.hpp:572
virtual returnValue setup()
virtual returnValue solveOutputSystem(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &index3, const ExportIndex &tmp_index, const ExportVariable &Ah, bool DERIVATIVES=false)
virtual returnValue sensitivitiesOutputSystem(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &index3, const ExportIndex &index4, const ExportIndex &tmp_index1, const ExportIndex &tmp_index2, const ExportVariable &Ah, const ExportVariable &Bh, bool STATES, uint number)
ExportVariable getRow(const ExportIndex &idx) const
ExportVariable rk_diffK
Definition: irk_export.hpp:571
double getFirstTime() const
DiagonallyImplicitRKExport(UserInteraction *_userInteraction=0, const std::string &_commonHeaderName="")
Definition: dirk_export.cpp:48
virtual returnValue evaluateMatrix(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &tmp_index, const ExportVariable &_rk_A, const ExportVariable &Ah, const ExportVariable &C, bool evaluateB, bool DERIVATIVES)
virtual ExportVariable getGlobalExportVariable(const uint factor) const
Allows to pass back messages to the calling function.
ExportVariable rk_kkk
Definition: rk_export.hpp:189
ExportLinearSolver * solver
Definition: irk_export.hpp:530
Block< Derived > block(Index startRow, Index startCol, Index blockRows, Index blockCols)
Definition: BlockMethods.h:56
BEGIN_NAMESPACE_ACADO typedef unsigned int uint
Definition: acado_types.hpp:42
virtual returnValue setup()
Allows to export code of a for-loop.
string toString(T const &value)
ExportVariable rk_mat1
Definition: irk_export.hpp:558
const std::string getNameOutputDiffs() const
ExportVariable getElement(const ExportIndex &rowIdx, const ExportIndex &colIdx) const
ExportVariable rk_diffsTemp2
ForwardIRKExport & operator=(const ForwardIRKExport &arg)
#define CLOSE_NAMESPACE_ACADO
ExportVariable rk_diffsNew3
DiagonallyImplicitRKExport & operator=(const DiagonallyImplicitRKExport &arg)
Definition: dirk_export.cpp:71
ExportVariable getSubMatrix(const ExportIndex &rowIdx1, const ExportIndex &rowIdx2, const ExportIndex &colIdx1, const ExportIndex &colIdx2) const
const std::string getNameDiffsRHS() const
Defines a scalar-valued index variable to be used for exporting code.
virtual returnValue evaluateStatesImplicitSystem(ExportStatementBlock *block, const ExportVariable &Ah, const ExportVariable &C, const ExportIndex &stage, const ExportIndex &i, const ExportIndex &j)
virtual ~DiagonallyImplicitRKExport()
Definition: dirk_export.cpp:61
ExportVariable rk_eta
Expression jacobian(const Expression &arg1, const Expression &arg2)
ExportStruct
static std::string fcnPrefix
DVector evaluateDerivedPolynomial(double time)
virtual DMatrix formMatrix(const DMatrix &mass, const DMatrix &jacobian)
virtual returnValue solveInputSystem(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &index3, const ExportIndex &tmp_index, const ExportVariable &Ah)
Definition: dirk_export.cpp:81
const std::string getNameOutputRHS() const
ExportVariable getCols(const ExportIndex &idx1, const ExportIndex &idx2) const
ExportVariable rk_auxSolver
Definition: irk_export.hpp:563
int acadoRoundAway(double x)
const std::string getNameSolveReuseFunction()
const std::string getNameRHS() const
const std::string get(const ExportIndex &rowIdx, const ExportIndex &colIdx) const
ExportVariable reset_int
virtual returnValue solveImplicitSystem(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &index3, const ExportIndex &tmp_index, const ExportVariable &Ah, const ExportVariable &C, const ExportVariable &det, bool DERIVATIVES=false)
Encapsulates all user interaction for setting options, logging data and plotting results.
virtual uint getDim() const
returnValue addStatement(const ExportStatement &_statement)
std::string getFullName() const
returnValue evaluateStatesOutputSystem(ExportStatementBlock *block, const ExportVariable &Ah, const ExportIndex &stage)
Definition: irk_export.cpp:705
unsigned getNumRows() const
Definition: matrix.hpp:185
ExportVariable rk_mat3
Definition: irk_export.hpp:567
ExportVariable rk_diffsNew2
unsigned getNumCols() const
Definition: matrix.hpp:189
uint getNumIntervals() const
virtual returnValue setup()=0
virtual returnValue evaluateRhsImplicitSystem(ExportStatementBlock *block, const ExportIndex &stage)
ExportVariable rk_xxx
std::string getName() const
ExportVariable getRows(const ExportIndex &idx1, const ExportIndex &idx2) const
returnValue addDeclaration(const ExportVariable &_data, ExportStruct _dataStruct=ACADO_ANY)
ExportAcadoFunction lin_input
double getLastTime() const
const std::string getNameSolveFunction()
#define BEGIN_NAMESPACE_ACADO
ExportVariable rk_rhsTemp
Definition: irk_export.hpp:564
virtual returnValue prepareInputSystem(ExportStatementBlock &code)
Allows to export a tailored implicit Runge-Kutta integrator with forward sensitivity generation for f...
ExportVariable rk_diffsTemp3
Definition: irk_export.hpp:569
virtual returnValue clear()
virtual returnValue prepareOutputSystem(ExportStatementBlock &code)
Allows to export code for a block of statements.
virtual returnValue sensitivitiesImplicitSystem(ExportStatementBlock *block, const ExportIndex &index1, const ExportIndex &index2, const ExportIndex &index3, const ExportIndex &tmp_index1, const ExportIndex &tmp_index2, const ExportVariable &Ah, const ExportVariable &Bh, const ExportVariable &det, bool STATES, uint number)
ExportArgument getAddress(const ExportIndex &_rowIdx, const ExportIndex &_colIdx=emptyConstExportIndex) const
GenericMatrix getRows(unsigned _start, unsigned _end) const
Definition: matrix.hpp:235
returnValue init(const uint newDim, const bool &reuse=true, const bool &unrolling=false)
ExportVariable getCol(const ExportIndex &idx) const
GenericVector< T > getRow(unsigned _idx) const
Definition: matrix.hpp:197
#define ACADOERROR(retval)
virtual bool equidistantControlGrid() const
Defines a matrix-valued variable to be used for exporting code.
returnValue addFunctionCall(const std::string &_fName, const ExportArgument &_argument1=emptyConstExportArgument, const ExportArgument &_argument2=emptyConstExportArgument, const ExportArgument &_argument3=emptyConstExportArgument, const ExportArgument &_argument4=emptyConstExportArgument, const ExportArgument &_argument5=emptyConstExportArgument, const ExportArgument &_argument6=emptyConstExportArgument, const ExportArgument &_argument7=emptyConstExportArgument, const ExportArgument &_argument8=emptyConstExportArgument, const ExportArgument &_argument9=emptyConstExportArgument)
Allows to export a tailored diagonally implicit Runge-Kutta integrator for fast model predictive cont...
Definition: dirk_export.hpp:54
std::string getName() const
Definition: export_data.cpp:86


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:32