source: src/main/java/agents/anac/y2019/harddealer/math3/ode/MultistepIntegrator.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 19.0 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package agents.anac.y2019.harddealer.math3.ode;
19
20import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
21import agents.anac.y2019.harddealer.math3.exception.MathIllegalStateException;
22import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
23import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
24import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
25import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
26import agents.anac.y2019.harddealer.math3.linear.Array2DRowRealMatrix;
27import agents.anac.y2019.harddealer.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
28import agents.anac.y2019.harddealer.math3.ode.nonstiff.DormandPrince853Integrator;
29import agents.anac.y2019.harddealer.math3.ode.sampling.StepHandler;
30import agents.anac.y2019.harddealer.math3.ode.sampling.StepInterpolator;
31import agents.anac.y2019.harddealer.math3.util.FastMath;
32
33/**
34 * This class is the base class for multistep integrators for Ordinary
35 * Differential Equations.
36 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
37 * <pre>
38 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
39 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
40 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
41 * ...
42 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
43 * </pre></p>
44 * <p>Rather than storing several previous steps separately, this implementation uses
45 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
46 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
47 * <pre>
48 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
49 * </pre>
50 * (we omit the k index in the notation for clarity)</p>
51 * <p>
52 * Multistep integrators with Nordsieck representation are highly sensitive to
53 * large step changes because when the step is multiplied by factor a, the
54 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
55 * and the last components are the least accurate ones. The default max growth
56 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
57 * </p>
58 *
59 * @see agents.anac.y2019.harddealer.math3.ode.nonstiff.AdamsBashforthIntegrator
60 * @see agents.anac.y2019.harddealer.math3.ode.nonstiff.AdamsMoultonIntegrator
61 * @since 2.0
62 */
63public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
64
65 /** First scaled derivative (h y'). */
66 protected double[] scaled;
67
68 /** Nordsieck matrix of the higher scaled derivatives.
69 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
70 */
71 protected Array2DRowRealMatrix nordsieck;
72
73 /** Starter integrator. */
74 private FirstOrderIntegrator starter;
75
76 /** Number of steps of the multistep method (excluding the one being computed). */
77 private final int nSteps;
78
79 /** Stepsize control exponent. */
80 private double exp;
81
82 /** Safety factor for stepsize control. */
83 private double safety;
84
85 /** Minimal reduction factor for stepsize control. */
86 private double minReduction;
87
88 /** Maximal growth factor for stepsize control. */
89 private double maxGrowth;
90
91 /**
92 * Build a multistep integrator with the given stepsize bounds.
93 * <p>The default starter integrator is set to the {@link
94 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
95 * some defaults settings.</p>
96 * <p>
97 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
98 * </p>
99 * @param name name of the method
100 * @param nSteps number of steps of the multistep method
101 * (excluding the one being computed)
102 * @param order order of the method
103 * @param minStep minimal step (must be positive even for backward
104 * integration), the last step can be smaller than this
105 * @param maxStep maximal step (must be positive even for backward
106 * integration)
107 * @param scalAbsoluteTolerance allowed absolute error
108 * @param scalRelativeTolerance allowed relative error
109 * @exception NumberIsTooSmallException if number of steps is smaller than 2
110 */
111 protected MultistepIntegrator(final String name, final int nSteps,
112 final int order,
113 final double minStep, final double maxStep,
114 final double scalAbsoluteTolerance,
115 final double scalRelativeTolerance)
116 throws NumberIsTooSmallException {
117
118 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
119
120 if (nSteps < 2) {
121 throw new NumberIsTooSmallException(
122 LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
123 nSteps, 2, true);
124 }
125
126 starter = new DormandPrince853Integrator(minStep, maxStep,
127 scalAbsoluteTolerance,
128 scalRelativeTolerance);
129 this.nSteps = nSteps;
130
131 exp = -1.0 / order;
132
133 // set the default values of the algorithm control parameters
134 setSafety(0.9);
135 setMinReduction(0.2);
136 setMaxGrowth(FastMath.pow(2.0, -exp));
137
138 }
139
140 /**
141 * Build a multistep integrator with the given stepsize bounds.
142 * <p>The default starter integrator is set to the {@link
143 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
144 * some defaults settings.</p>
145 * <p>
146 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
147 * </p>
148 * @param name name of the method
149 * @param nSteps number of steps of the multistep method
150 * (excluding the one being computed)
151 * @param order order of the method
152 * @param minStep minimal step (must be positive even for backward
153 * integration), the last step can be smaller than this
154 * @param maxStep maximal step (must be positive even for backward
155 * integration)
156 * @param vecAbsoluteTolerance allowed absolute error
157 * @param vecRelativeTolerance allowed relative error
158 */
159 protected MultistepIntegrator(final String name, final int nSteps,
160 final int order,
161 final double minStep, final double maxStep,
162 final double[] vecAbsoluteTolerance,
163 final double[] vecRelativeTolerance) {
164 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
165 starter = new DormandPrince853Integrator(minStep, maxStep,
166 vecAbsoluteTolerance,
167 vecRelativeTolerance);
168 this.nSteps = nSteps;
169
170 exp = -1.0 / order;
171
172 // set the default values of the algorithm control parameters
173 setSafety(0.9);
174 setMinReduction(0.2);
175 setMaxGrowth(FastMath.pow(2.0, -exp));
176
177 }
178
179 /**
180 * Get the starter integrator.
181 * @return starter integrator
182 */
183 public ODEIntegrator getStarterIntegrator() {
184 return starter;
185 }
186
187 /**
188 * Set the starter integrator.
189 * <p>The various step and event handlers for this starter integrator
190 * will be managed automatically by the multi-step integrator. Any
191 * user configuration for these elements will be cleared before use.</p>
192 * @param starterIntegrator starter integrator
193 */
194 public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
195 this.starter = starterIntegrator;
196 }
197
198 /** Start the integration.
199 * <p>This method computes one step using the underlying starter integrator,
200 * and initializes the Nordsieck vector at step start. The starter integrator
201 * purpose is only to establish initial conditions, it does not really change
202 * time by itself. The top level multistep integrator remains in charge of
203 * handling time propagation and events handling as it will starts its own
204 * computation right from the beginning. In a sense, the starter integrator
205 * can be seen as a dummy one and so it will never trigger any user event nor
206 * call any user step handler.</p>
207 * @param t0 initial time
208 * @param y0 initial value of the state vector at t0
209 * @param t target time for the integration
210 * (can be set to a value smaller than <code>t0</code> for backward integration)
211 * @exception DimensionMismatchException if arrays dimension do not match equations settings
212 * @exception NumberIsTooSmallException if integration step is too small
213 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
214 * @exception NoBracketingException if the location of an event cannot be bracketed
215 */
216 protected void start(final double t0, final double[] y0, final double t)
217 throws DimensionMismatchException, NumberIsTooSmallException,
218 MaxCountExceededException, NoBracketingException {
219
220 // make sure NO user event nor user step handler is triggered,
221 // this is the task of the top level integrator, not the task
222 // of the starter integrator
223 starter.clearEventHandlers();
224 starter.clearStepHandlers();
225
226 // set up one specific step handler to extract initial Nordsieck vector
227 starter.addStepHandler(new NordsieckInitializer((nSteps + 3) / 2, y0.length));
228
229 // start integration, expecting a InitializationCompletedMarkerException
230 try {
231
232 if (starter instanceof AbstractIntegrator) {
233 ((AbstractIntegrator) starter).integrate(getExpandable(), t);
234 } else {
235 starter.integrate(new FirstOrderDifferentialEquations() {
236
237 /** {@inheritDoc} */
238 public int getDimension() {
239 return getExpandable().getTotalDimension();
240 }
241
242 /** {@inheritDoc} */
243 public void computeDerivatives(double t, double[] y, double[] yDot) {
244 getExpandable().computeDerivatives(t, y, yDot);
245 }
246
247 }, t0, y0, t, new double[y0.length]);
248 }
249
250 // we should not reach this step
251 throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
252
253 } catch (InitializationCompletedMarkerException icme) { // NOPMD
254 // this is the expected nominal interruption of the start integrator
255
256 // count the evaluations used by the starter
257 getCounter().increment(starter.getEvaluations());
258
259 }
260
261 // remove the specific step handler
262 starter.clearStepHandlers();
263
264 }
265
266 /** Initialize the high order scaled derivatives at step start.
267 * @param h step size to use for scaling
268 * @param t first steps times
269 * @param y first steps states
270 * @param yDot first steps derivatives
271 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
272 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
273 */
274 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
275 final double[][] y,
276 final double[][] yDot);
277
278 /** Get the minimal reduction factor for stepsize control.
279 * @return minimal reduction factor
280 */
281 public double getMinReduction() {
282 return minReduction;
283 }
284
285 /** Set the minimal reduction factor for stepsize control.
286 * @param minReduction minimal reduction factor
287 */
288 public void setMinReduction(final double minReduction) {
289 this.minReduction = minReduction;
290 }
291
292 /** Get the maximal growth factor for stepsize control.
293 * @return maximal growth factor
294 */
295 public double getMaxGrowth() {
296 return maxGrowth;
297 }
298
299 /** Set the maximal growth factor for stepsize control.
300 * @param maxGrowth maximal growth factor
301 */
302 public void setMaxGrowth(final double maxGrowth) {
303 this.maxGrowth = maxGrowth;
304 }
305
306 /** Get the safety factor for stepsize control.
307 * @return safety factor
308 */
309 public double getSafety() {
310 return safety;
311 }
312
313 /** Set the safety factor for stepsize control.
314 * @param safety safety factor
315 */
316 public void setSafety(final double safety) {
317 this.safety = safety;
318 }
319
320 /** Get the number of steps of the multistep method (excluding the one being computed).
321 * @return number of steps of the multistep method (excluding the one being computed)
322 */
323 public int getNSteps() {
324 return nSteps;
325 }
326
327 /** Compute step grow/shrink factor according to normalized error.
328 * @param error normalized error of the current step
329 * @return grow/shrink factor for next step
330 */
331 protected double computeStepGrowShrinkFactor(final double error) {
332 return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
333 }
334
335 /** Transformer used to convert the first step to Nordsieck representation.
336 * @deprecated as of 3.6 this unused interface is deprecated
337 */
338 @Deprecated
339 public interface NordsieckTransformer {
340 /** Initialize the high order scaled derivatives at step start.
341 * @param h step size to use for scaling
342 * @param t first steps times
343 * @param y first steps states
344 * @param yDot first steps derivatives
345 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
346 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
347 */
348 Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
349 final double[][] y,
350 final double[][] yDot);
351 }
352
353 /** Specialized step handler storing the first step. */
354 private class NordsieckInitializer implements StepHandler {
355
356 /** Steps counter. */
357 private int count;
358
359 /** First steps times. */
360 private final double[] t;
361
362 /** First steps states. */
363 private final double[][] y;
364
365 /** First steps derivatives. */
366 private final double[][] yDot;
367
368 /** Simple constructor.
369 * @param nbStartPoints number of start points (including the initial point)
370 * @param n problem dimension
371 */
372 NordsieckInitializer(final int nbStartPoints, final int n) {
373 this.count = 0;
374 this.t = new double[nbStartPoints];
375 this.y = new double[nbStartPoints][n];
376 this.yDot = new double[nbStartPoints][n];
377 }
378
379 /** {@inheritDoc} */
380 public void handleStep(StepInterpolator interpolator, boolean isLast)
381 throws MaxCountExceededException {
382
383 final double prev = interpolator.getPreviousTime();
384 final double curr = interpolator.getCurrentTime();
385
386 if (count == 0) {
387 // first step, we need to store also the point at the beginning of the step
388 interpolator.setInterpolatedTime(prev);
389 t[0] = prev;
390 final ExpandableStatefulODE expandable = getExpandable();
391 final EquationsMapper primary = expandable.getPrimaryMapper();
392 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
393 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
394 int index = 0;
395 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
396 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
397 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
398 ++index;
399 }
400 }
401
402 // store the point at the end of the step
403 ++count;
404 interpolator.setInterpolatedTime(curr);
405 t[count] = curr;
406
407 final ExpandableStatefulODE expandable = getExpandable();
408 final EquationsMapper primary = expandable.getPrimaryMapper();
409 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
410 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
411 int index = 0;
412 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
413 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
414 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
415 ++index;
416 }
417
418 if (count == t.length - 1) {
419
420 // this was the last point we needed, we can compute the derivatives
421 stepStart = t[0];
422 stepSize = (t[t.length - 1] - t[0]) / (t.length - 1);
423
424 // first scaled derivative
425 scaled = yDot[0].clone();
426 for (int j = 0; j < scaled.length; ++j) {
427 scaled[j] *= stepSize;
428 }
429
430 // higher order derivatives
431 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
432
433 // stop the integrator now that all needed steps have been handled
434 throw new InitializationCompletedMarkerException();
435
436 }
437
438 }
439
440 /** {@inheritDoc} */
441 public void init(double t0, double[] y0, double time) {
442 // nothing to do
443 }
444
445 }
446
447 /** Marker exception used ONLY to stop the starter integrator after first step. */
448 private static class InitializationCompletedMarkerException
449 extends RuntimeException {
450
451 /** Serializable version identifier. */
452 private static final long serialVersionUID = -1914085471038046418L;
453
454 /** Simple constructor. */
455 InitializationCompletedMarkerException() {
456 super((Throwable) null);
457 }
458
459 }
460
461}
Note: See TracBrowser for help on using the repository browser.