source: branches/dev-5819/Python/src/tests_jmodelica/optimization/test_casadi_collocation.py @ 13600

Last change on this file since 13600 was 13600, checked in by randersson, 2 months ago

#5819 Merged JM trunk into branch.

File size: 90.1 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4# Copyright (C) 2014 Modelon AB
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation, version 3 of the License.
9#
10# This program is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13# GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License
16# along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18"""Tests the casadi_collocation module."""
19
20import os
21import nose
22
23from collections import OrderedDict
24import numpy as N
25from scipy.io.matlab.mio import loadmat
26
27from tests_jmodelica import testattr, get_files_path
28from pyjmi.common.io import ResultDymolaTextual
29from pymodelica import compile_fmu
30from pyfmi import load_fmu
31try:
32    from pyjmi import transfer_to_casadi_interface
33    from pyjmi.optimization.casadi_collocation import *
34    import casadi
35except (NameError, ImportError):
36    pass
37
38from pyjmi.common.io import VariableNotFoundError as jmiVariableNotFoundError
39#Check to see if pyfmi is installed so that we also catch the error generated
40#from that package
41try:
42    from pyfmi.common.io import VariableNotFoundError as fmiVariableNotFoundError
43    VariableNotFoundError = (jmiVariableNotFoundError, fmiVariableNotFoundError)
44except ImportError:
45    VariableNotFoundError = jmiVariableNotFoundError
46
47path_to_mos = os.path.join(get_files_path(), 'Modelica')
48path_to_data = os.path.join(get_files_path(), 'Data')
49
50def assert_results(res, cost_ref, u_norm_ref,
51                   cost_rtol=1e-3, u_norm_rtol=1e-4, input_name="u"):
52    """Helper function for asserting optimization results."""
53    cost = float(res.solver.solver_object.output(casadi.NLP_SOLVER_F))
54    u = res[input_name]
55    u_norm = N.linalg.norm(u) / N.sqrt(len(u))
56    N.testing.assert_allclose(cost, cost_ref, cost_rtol)
57    N.testing.assert_allclose(u_norm, u_norm_ref, u_norm_rtol)
58
59class TestLocalDAECollocator(object):
60   
61    """
62    Tests pyjmi.optimization.casadi_collocation.LocalDAECollocator.
63    """
64   
65    @classmethod
66    def setUpClass(self):
67        """Compile the test models."""
68        self.compiler_opts = self.compile_options()
69        vdp_file_path = os.path.join(get_files_path(), 'Modelica', 'VDP.mop')
70        class_path = "VDP_pack.VDP_Opt_Bounds_Lagrange"
71        self.vdp_bounds_lagrange_op = \
72                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
73       
74        class_path = "VDP_pack.VDP_Opt_Bounds_Lagrange_Renamed_Input"
75        self.vdp_bounds_lagrange_renamed_op = \
76                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
77       
78        class_path = "VDP_pack.VDP_Opt_Bounds_Mayer"
79        self.vdp_bounds_mayer_op = \
80                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
81       
82        class_path = "VDP_pack.VDP_Opt_Constraints_Mayer"
83        self.vdp_constraints_mayer_op = \
84                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
85       
86        class_path = "VDP_pack.VDP_Opt_Initial_Equations"
87        self.vdp_initial_equations_op = \
88                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
89       
90        class_path = "VDP_pack.VDP_Opt_Scaled_Min_Time"
91        self.vdp_scaled_min_time_op = \
92                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
93       
94        class_path = "VDP_pack.VDP_Opt_Unscaled_Min_Time"
95        self.vdp_unscaled_min_time_op = \
96                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
97
98        class_path = "VDP_pack.VDP_Opt_Min_Time_Nonzero_Start"
99        self.vdp_min_time_nonzero_start_op = \
100                transfer_to_casadi_interface(class_path, vdp_file_path, self.compiler_opts)
101       
102        class_path = "VDP_pack.VDP_Opt_Function"
103        self.vdp_function_op = \
104                transfer_to_casadi_interface(
105                        class_path, vdp_file_path,
106                        compiler_options=dict(list({"inline_functions": "none"}.items()) + list(self.compiler_opts.items())))
107       
108        class_path = "VDP_pack.VDP_Opt_Function_global_constant"
109        self.vdp_function_op_global_constant = \
110                transfer_to_casadi_interface(
111                        class_path, vdp_file_path,
112                        compiler_options=dict(list({"inline_functions": "none"}.items()) + list(self.compiler_opts.items())))
113       
114        cstr_file_path = os.path.join(get_files_path(), 'Modelica', 'CSTR.mop')
115        class_path = "CSTR.CSTR"
116        fmu_cstr = compile_fmu(class_path, cstr_file_path,
117                               separate_process=True)
118        self.cstr_model = load_fmu(fmu_cstr)
119       
120        class_path = "CSTR.CSTR_Opt_Bounds_Lagrange"
121        self.cstr_lagrange_op = \
122                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
123       
124        class_path = "CSTR.CSTR_Opt_Bounds_Mayer"
125        self.cstr_mayer_op = \
126                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
127       
128        class_path = "CSTR.CSTR_Opt_Dependent_Parameter"
129        self.cstr_dependent_parameter_op = \
130                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
131       
132        class_path = "CSTR.CSTR_Opt_Extends"
133        self.cstr_extends_op = \
134                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
135       
136        class_path = "CSTR.CSTR_Opt_Scaled_Min_Time"
137        self.cstr_scaled_min_time_op = \
138                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
139       
140        class_path = "CSTR.CSTR_Opt_Unscaled_Min_Time"
141        self.cstr_unscaled_min_time_op = \
142                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
143
144        class_path = "CSTR.CSTR_Opt_Min_Time_No_Initial_Guess"
145        self.cstr_min_time_no_init_guess_op = \
146                transfer_to_casadi_interface(class_path, cstr_file_path, self.compiler_opts)
147       
148        pe_file_path = os.path.join(get_files_path(), 'Modelica',
149                                    'ParameterEstimation_1.mop')
150        class_path = "ParEst.SecondOrder"
151        fmu_second_order = compile_fmu(class_path, pe_file_path,
152                                       separate_process=True)
153        self.second_order_model = load_fmu(fmu_second_order)
154       
155        class_path = "ParEst.ParEstCasADi"
156        self.second_order_par_est_op = \
157                transfer_to_casadi_interface(class_path, pe_file_path, self.compiler_opts)
158       
159        qt_file_path = os.path.join(get_files_path(), 'Modelica',
160                                    'QuadTankPack.mop')
161        class_path = "QuadTankPack.Sim_QuadTank"
162        fmu_qt_sim = compile_fmu(class_path, qt_file_path,
163                                 separate_process=True)
164        self.qt_model = load_fmu(fmu_qt_sim)
165       
166        class_path = "QuadTankPack.QuadTank_ParEstCasADi"
167        self.qt_par_est_op = \
168                transfer_to_casadi_interface(class_path, qt_file_path, self.compiler_opts)
169       
170        class_path = "QuadTankPack.QuadTank_ParEstCasADi_Degenerate"
171        self.qt_par_est_degenerate_op = \
172                transfer_to_casadi_interface(class_path, qt_file_path, self.compiler_opts)
173               
174        stream_file_path = os.path.join(get_files_path(), 'Modelica', 'TestOptimizationProblems.mop')
175        class_path = "InStream"
176        stream_comp_opts = self.compiler_opts
177        stream_comp_opts['eliminate_alias_variables'] = False
178        stream_comp_opts['eliminate_linear_equations'] = False
179        stream_comp_opts['variability_propagation'] = False
180        self.stream_op = transfer_to_casadi_interface(class_path, stream_file_path, stream_comp_opts)
181       
182        self.algorithm = "LocalDAECollocationAlg"
183       
184    def optimize_options(self, op, alg=None):
185        if alg is None:
186            return op.optimize_options()
187        else:
188            return op.optimize_options(alg)
189   
190    @staticmethod
191    def compile_options():
192        compiler_options = {"common_subexp_elim":False}
193        return compiler_options
194
195    @testattr(casadi_base = True)
196    def test_free_horizon_tdlagrange(self):
197        """
198        Test solving a minimum time problem with a time dependent Lagrange
199        integrand.
200        """
201       
202        op = self.vdp_unscaled_min_time_op
203        objInt = op.getObjectiveIntegrand()
204        op.setObjectiveIntegrand(op.getTimeVariable())
205       
206        # References values
207        tf = 2.2811590707107996e0  # old cost_ref
208        cost_ref =  tf + (tf**2)/2
209        u_norm_ref = 9.991517452037317e-1
210
211        # Scaled, Radau
212        opts = self.optimize_options(op, self.algorithm)
213        opts['discr'] = "LGR"
214        res = op.optimize(self.algorithm, opts)
215        assert_results(res, cost_ref, u_norm_ref)
216       
217        # With checkpoints
218        opts['checkpoint'] = True
219        opts['expand_to_sx'] = "DAE"
220        res = op.optimize(self.algorithm, opts)
221        assert_results(res, cost_ref, u_norm_ref)
222        op.setObjectiveIntegrand(objInt)
223       
224    @testattr(casadi_base = True)
225    def test_result_modes(self):
226        """
227        Test all combinations of discretization schemes and result modes.
228        """
229        cost_ref = 1.8576873858261e3
230        u_norm_ref = 3.050971000653911e2
231       
232        op = self.cstr_lagrange_op
233       
234        # Shift time horizon
235        t0 = op.get('startTime')
236        tf = op.get('finalTime')
237        dt = 5.
238        op.set('startTime', t0 + dt)
239        op.set('finalTime', tf + dt)
240
241        # Define function for checking time horizon
242        def assert_time(res):
243            N.testing.assert_allclose(res['time'][[0, -1]], [t0 + dt, tf + dt])
244
245        opts = self.optimize_options(op, self.algorithm)
246       
247        opts['discr'] = 'LG'
248        opts['result_mode'] = 'collocation_points'
249        print('LG + collocation_points')
250        res = op.optimize(self.algorithm, opts)
251        assert_time(res)
252        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
253       
254        opts['discr'] = 'LG'
255        opts['result_mode'] = 'element_interpolation'
256        print('LG + element_interpolation')
257        res = op.optimize(self.algorithm, opts)
258        N.testing.assert_allclose(res['time'][[0, -1]], [t0 + dt, tf + dt])
259        assert_time(res)
260        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
261       
262        opts['discr'] = 'LG'
263        opts['result_mode'] = 'mesh_points'
264        print('LG + mesh_points')
265        res = op.optimize(self.algorithm, opts)
266        assert_time(res)
267        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
268       
269        opts['discr'] = 'LGR'
270        opts['result_mode'] = 'collocation_points'
271        print('LGR + collocation_points')
272        res = op.optimize(self.algorithm, opts)
273        assert_time(res)
274        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
275       
276        opts['discr'] = 'LGR'
277        opts['result_mode'] = 'element_interpolation'
278        print('LGR + element_interpolation')
279        res = op.optimize(self.algorithm, opts)
280        assert_time(res)
281        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
282       
283        opts['discr'] = 'LGR'
284        opts['result_mode'] = 'mesh_points'
285        print('LGR + mesh_points')
286        res = op.optimize(self.algorithm, opts)
287        assert_time(res)
288        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
289
290        # Reset time horizon
291        op.set('startTime', t0 - dt)
292        op.set('finalTime', tf - dt)
293       
294    @testattr(casadi_base = True)
295    def test_nominal_zero(self):
296        """
297        Test a parameter estimation example with a nominal trajectory equal to
298        zero and the same example with a nominal attribute equal to zero.
299        """
300        model = self.second_order_model
301        model.reset()
302        op = self.second_order_par_est_op
303       
304        # Simulate with initial guesses
305        model.set('w', 1.3)
306        model.set('z', 0)
307        sim_res = model.simulate(final_time=15.)
308       
309        # Reference values
310        w_ref = 1.048589
311        z_ref = 0.470934
312       
313        # Measurements 
314        y_meas = N.array([0.01463904, 0.35424225, 0.94776816, 1.20116167,
315                          1.17283905, 1.03631145, 1.0549561, 0.94827652,
316                          1.0317119, 1.04010453, 1.08012155])
317        t_meas = N.linspace(0., 10., num=len(y_meas))
318        data = N.vstack([t_meas, y_meas])
319       
320        # Measurement data
321        Q = N.array([[1.]])
322        quad_pen = OrderedDict()
323        quad_pen['y'] = data
324        external_data = ExternalData(quad_pen=quad_pen, Q=Q)
325       
326        # Optimization options
327        opts = self.optimize_options(op, self.algorithm)
328        opts['external_data'] = external_data
329        opts['n_e'] = 16
330       
331        # Optimize with trajectories
332        opts['nominal_traj'] = sim_res
333        traj_res = op.optimize(self.algorithm, opts)
334        w_traj = traj_res.final('w')
335        z_traj = traj_res.final('z')
336        N.testing.assert_allclose(w_traj, w_ref, 1e-2)
337        N.testing.assert_allclose(z_traj, z_ref, 1e-2)
338       
339        # Set nominal attribute to zero
340        op.getVariable('z').setNominal(0)
341
342        # Optimize with scaling
343        opts = self.optimize_options(op, self.algorithm)
344        opts['external_data'] = external_data
345        opts['variable_scaling'] = True
346        N.testing.assert_raises(CasadiCollocatorException,
347                                op.optimize, self.algorithm, opts)
348        op.getVariable('z').setNominal(1)
349       
350    @testattr(casadi_base = True)
351    def test_init_traj_sim(self):
352        """Test initial trajectories based on an existing simulation."""
353        model = self.cstr_model
354        model.reset()
355        op = self.cstr_extends_op
356        op_min_time = self.cstr_min_time_no_init_guess_op
357        model.set(['c_init', 'T_init'], op.get(['c_init', 'T_init']))
358       
359        # Create input trajectory
360        t = [0, 200]
361        u = [342.85, 280]
362        u_traj = N.transpose(N.vstack((t, u)))
363       
364        # Generate initial trajectories
365        opts = model.simulate_options()
366        opts["CVode_options"]["rtol"] = 1e-6
367        opts["CVode_options"]["atol"] = 1e-8 * model.nominal_continuous_states
368        init_res = model.simulate(final_time=300, input=('Tc', u_traj),
369                                  options=opts)
370       
371        # Check computed initial guess
372        opts = self.optimize_options(op, self.algorithm)
373        opts['variable_scaling'] = False
374        opts['init_traj'] = init_res
375        col = LocalDAECollocator(op, opts)
376        xx_init = col.get_xx_init()
377        N.testing.assert_allclose(
378                xx_init[col.var_indices['x'][opts['n_e']][opts['n_cp']]],
379                [435.4425832, 333.42862629], rtol=1e-4)
380
381        # Check compute initial guess for minimum time
382        opts = self.optimize_options(op_min_time, self.algorithm)
383        opts['variable_scaling'] = False
384        opts['init_traj'] = init_res
385        col = LocalDAECollocator(op_min_time, opts)
386        xx_init = col.get_xx_init()
387        (ind, vt) = col.name_map['finalTime']
388        N.testing.assert_allclose(xx_init[col.var_indices[vt][ind]], 300.,
389                                  rtol=1e-4)
390
391    @testattr(casadi_base = True)
392    def test_init_dual(self):
393        """Test initializing dual variables."""
394        op = self.vdp_bounds_lagrange_op
395       
396        # References values
397        cost_ref = 3.19495079586595e0
398        u_norm_ref = 2.80997269112246e-1
399       
400        # Get initial guess
401        opts = self.optimize_options(op, self.algorithm)
402        opts['n_e'] = 40
403        opts['n_cp'] = 2
404        res = op.optimize(self.algorithm, opts)
405        assert_results(res, cost_ref, u_norm_ref)
406       
407        # Optimize without warm start
408        opts['init_traj'] = res
409        opts['init_dual'] = res.dual_opt
410        res = op.optimize(self.algorithm, opts)
411        assert_results(res, cost_ref, u_norm_ref, 5e-2, 5e-2)
412
413        # Optimize with warm start
414        opts['IPOPT_options']['warm_start_init_point'] = "yes"
415        opts['IPOPT_options']['mu_init'] = 10 ** (-8.6)
416        opts['IPOPT_options']['warm_start_bound_push'] = 1e-12
417        opts['IPOPT_options']['warm_start_bound_frac'] = 1e-12
418        opts['IPOPT_options']['warm_start_slack_bound_frac'] = 1e-12
419        opts['IPOPT_options']['warm_start_slack_bound_push'] = 1e-12
420        opts['IPOPT_options']['warm_start_mult_bound_push'] = 1e-12
421        res = op.optimize(self.algorithm, opts)
422        assert_results(res, cost_ref, u_norm_ref, 5e-2, 5e-2)
423
424    @testattr(casadi_base = True)
425    def test_init_traj_opt(self):
426        """Test optimizing based on an existing optimization reult."""
427        op = self.vdp_bounds_lagrange_op
428       
429        # References values
430        cost_ref = 3.19495079586595e0
431        u_norm_ref = 2.80997269112246e-1
432       
433        # Get initial guess
434        opts = self.optimize_options(op, self.algorithm)
435        opts['n_e'] = 40
436        opts['n_cp'] = 2
437        res = op.optimize(self.algorithm, opts)
438        assert_results(res, cost_ref, u_norm_ref)
439       
440        # Optimize using initial guess
441        opts['n_e'] = 75
442        opts['n_cp'] = 4
443        #~ opts['eliminate_der_var'] = True #3509
444        #~ opts['eliminate_cont_var'] = True #3509
445        opts['init_traj'] = ResultDymolaTextual(
446                "VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt")
447        res = op.optimize(self.algorithm, opts)
448        assert_results(res, cost_ref, u_norm_ref, 5e-2, 5e-2)
449   
450    """
451    @testattr(casadi_base = True)
452    def test_nominal_traj_with_check_point_vdp(self):
453        op = self.vdp_bounds_lagrange_op
454       
455        # References values
456        cost_ref_traj = 3.19495079586595e0
457        u_norm_ref_traj = 2.80997269112246e-1
458        cost_ref = 3.1749908234182826e0
459        u_norm_ref = 2.848606420347583e-1
460       
461        # Get nominal and initial trajectories
462        opts = self.optimize_options(op, self.algorithm)
463        opts['n_e'] = 40
464        opts['n_cp'] = 2
465        opts["checkpoint"] = True
466        res = op.optimize(self.algorithm, opts)
467        assert_results(res, cost_ref_traj, u_norm_ref_traj)
468
469        init_data = ResultDymolaTextual(res.result_file)
470
471        # Optimize using only initial trajectories
472        opts['n_e'] = 75
473        opts['n_cp'] = 4
474        opts['init_traj'] = init_data
475        res = op.optimize(self.algorithm, opts)
476        assert_results(res, cost_ref, u_norm_ref)
477       
478        # Optimize using nominal and initial trajectories
479        opts['nominal_traj'] = init_data
480        opts['nominal_traj_mode'] = {'_default_mode': "affine"}
481        res = op.optimize(self.algorithm, opts)
482        assert_results(res, cost_ref, u_norm_ref)
483        col = res.solver
484        xx_init = col.get_xx_init()
485        N.testing.assert_allclose(
486                xx_init[col.var_indices['x'][opts['n_e']][opts['n_cp']]],
487                [0.85693481, 0.12910473])
488
489        # Test disabling scaling
490        opts['variable_scaling'] = False
491        res = op.optimize(self.algorithm, opts)
492        assert_results(res, cost_ref, u_norm_ref)
493    """
494   
495    @testattr(casadi_base = True)
496    def test_updated_nominal_traj_vdp(self):
497        """Test optimizing a VDP using nominal and initial trajectories."""
498        op = self.vdp_bounds_lagrange_op
499       
500        # References values
501        cost_ref_traj = 3.19495079586595e0
502        u_norm_ref_traj = 2.80997269112246e-1
503        cost_ref = 3.1749908234182826e0
504        u_norm_ref = 2.848606420347583e-1
505       
506        # Get nominal and initial trajectories
507        opts = self.optimize_options(op, self.algorithm)
508        opts['n_e'] = 40
509        opts['n_cp'] = 2
510        opts["variable_scaling_allow_update"] = True
511        res = op.optimize(self.algorithm, opts)
512        assert_results(res, cost_ref_traj, u_norm_ref_traj)
513        try:
514            os.remove("vdp_nom_traj_result.txt")
515        except OSError:
516            pass
517        os.rename("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt",
518                  "vdp_nom_traj_result.txt")
519       
520        solver = res.get_solver()
521        solver.set_nominal_traj(ResultDymolaTextual("vdp_nom_traj_result.txt"), {'_default_mode': "affine"})
522       
523        res_update = solver.optimize()
524       
525        cost_update = float(res_update.solver.solver_object.output(casadi.NLP_SOLVER_F))
526       
527        # Optimize using nominal and initial trajectories
528        opts['nominal_traj'] = ResultDymolaTextual("vdp_nom_traj_result.txt")
529        opts['nominal_traj_mode'] = {'_default_mode': "affine"}
530        res = op.optimize(self.algorithm, opts)
531       
532        cost = float(res.solver.solver_object.output(casadi.NLP_SOLVER_F))
533
534        N.testing.assert_equal(cost_update, cost)
535       
536    @testattr(casadi_base = True)
537    def test_no_updated_nominal_traj_vdp(self):
538        """Test optimizing a VDP using nominal and initial trajectories."""
539        op = self.vdp_bounds_lagrange_op
540       
541        # References values
542        cost_ref_traj = 3.19495079586595e0
543        u_norm_ref_traj = 2.80997269112246e-1
544        cost_ref = 3.1749908234182826e0
545        u_norm_ref = 2.848606420347583e-1
546       
547        # Get nominal and initial trajectories
548        opts = self.optimize_options(op, self.algorithm)
549        opts['n_e'] = 40
550        opts['n_cp'] = 2
551        opts["variable_scaling_allow_update"] = False
552        res = op.optimize(self.algorithm, opts)
553        assert_results(res, cost_ref_traj, u_norm_ref_traj)
554        try:
555            os.remove("vdp_nom_traj_result.txt")
556        except OSError:
557            pass
558        os.rename("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt",
559                  "vdp_nom_traj_result.txt")
560       
561        assert res.get_solver().collocator.pp.shape[0] == 2 #Two parameters
562       
563        solver = res.get_solver()
564       
565        N.testing.assert_raises(CasadiCollocatorException,
566                                solver.set_nominal_traj, res)
567
568   
569    @testattr(casadi_base = True)
570    def test_nominal_traj_vdp(self):
571        """Test optimizing a VDP using nominal and initial trajectories."""
572        op = self.vdp_bounds_lagrange_op
573       
574        # References values
575        cost_ref_traj = 3.19495079586595e0
576        u_norm_ref_traj = 2.80997269112246e-1
577        cost_ref = 3.1749908234182826e0
578        u_norm_ref = 2.848606420347583e-1
579       
580        # Get nominal and initial trajectories
581        opts = self.optimize_options(op, self.algorithm)
582        opts['n_e'] = 40
583        opts['n_cp'] = 2
584        res = op.optimize(self.algorithm, opts)
585        assert_results(res, cost_ref_traj, u_norm_ref_traj)
586        try:
587            os.remove("vdp_nom_traj_result.txt")
588        except OSError:
589            pass
590        os.rename("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt",
591                  "vdp_nom_traj_result.txt")
592       
593        # Optimize using only initial trajectories
594        opts['n_e'] = 75
595        opts['n_cp'] = 4
596        opts['init_traj'] = ResultDymolaTextual("vdp_nom_traj_result.txt")
597        res = op.optimize(self.algorithm, opts)
598        assert_results(res, cost_ref, u_norm_ref)
599       
600        # Optimize using nominal and initial trajectories
601        opts['nominal_traj'] = ResultDymolaTextual("vdp_nom_traj_result.txt")
602        opts['nominal_traj_mode'] = {'_default_mode': "affine"}
603        res = op.optimize(self.algorithm, opts)
604        assert_results(res, cost_ref, u_norm_ref)
605        col = res.solver
606        xx_init = col.get_xx_init()
607        N.testing.assert_allclose(
608                xx_init[col.var_indices['x'][opts['n_e']][opts['n_cp']]],
609                [0.85693481, 0.12910473])
610
611        # #3509
612        #~ # Test with eliminated continuity variables
613        #~ opts['eliminate_cont_var'] = True
614        #~ res = op.optimize(self.algorithm, opts)
615        #~ assert_results(res, cost_ref, u_norm_ref)
616        #~
617        #~ # Test with eliminated continuity and derivative variables
618        #~ opts['eliminate_der_var'] = True
619        #~ res = op.optimize(self.algorithm, opts)
620        #~ assert_results(res, cost_ref, u_norm_ref)
621       
622        # Test disabling scaling
623        opts['variable_scaling'] = False
624        res = op.optimize(self.algorithm, opts)
625        assert_results(res, cost_ref, u_norm_ref)
626       
627    @testattr(casadi_base = True)
628    def test_nominal_traj_vdp_lg(self):
629        """Test optimizing a VDP using nominal and initial trajectories."""
630        op = self.vdp_bounds_lagrange_op
631       
632        # References values
633        #cost_ref_traj = 3.19495079586595e0
634        cost_ref_traj = 3.171799772935013
635        #u_norm_ref_traj = 2.80997269112246e-1
636        u_norm_ref_traj = 2.8094758146374343e-1
637        cost_ref = 3.1749908234182826e0
638        #u_norm_ref = 2.848606420347583e-1
639        u_norm_ref = 2.844516216499257e-1
640       
641        # Get nominal and initial trajectories
642        opts = self.optimize_options(op, self.algorithm)
643        opts['n_e'] = 40
644        opts['n_cp'] = 2
645        opts['discr'] = "LG" #TEST LG DISCR
646        res = op.optimize(self.algorithm, opts)
647        assert_results(res, cost_ref_traj, u_norm_ref_traj)
648        try:
649            os.remove("vdp_nom_traj_result.txt")
650        except OSError:
651            pass
652        os.rename("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt",
653                  "vdp_nom_traj_result.txt")
654       
655        # Optimize using only initial trajectories
656        opts['n_e'] = 75
657        opts['n_cp'] = 4
658        opts['init_traj'] = ResultDymolaTextual("vdp_nom_traj_result.txt")
659        res = op.optimize(self.algorithm, opts)
660        assert_results(res, cost_ref, u_norm_ref)
661       
662        # Optimize using nominal and initial trajectories
663        opts['nominal_traj'] = ResultDymolaTextual("vdp_nom_traj_result.txt")
664        opts['nominal_traj_mode'] = {'_default_mode': "affine"}
665        res = op.optimize(self.algorithm, opts)
666        assert_results(res, cost_ref, u_norm_ref)
667        col = res.solver
668        xx_init = col.get_xx_init()
669        N.testing.assert_allclose(
670                xx_init[col.var_indices['x'][opts['n_e']][opts['n_cp']]],
671                [0.866368, 0.118284], rtol=1e-5) #[0.85693481, 0.12910473]
672
673        # #3509
674        #~ # Test with eliminated continuity variables
675        #~ opts['eliminate_cont_var'] = True
676        #~ res = op.optimize(self.algorithm, opts)
677        #~ assert_results(res, cost_ref, u_norm_ref)
678        #~
679        #~ # Test with eliminated continuity and derivative variables
680        #~ opts['eliminate_der_var'] = True
681        #~ res = op.optimize(self.algorithm, opts)
682        #~ assert_results(res, cost_ref, u_norm_ref)
683       
684        # Test disabling scaling
685        opts['variable_scaling'] = False
686        res = op.optimize(self.algorithm, opts)
687        assert_results(res, cost_ref, u_norm_ref)
688
689    @testattr(casadi_base = True)
690    def test_nominal_traj_cstr(self):
691        """Test optimizing a CSTR using nominal and initial trajectories."""
692        op = self.cstr_lagrange_op
693       
694        # References values
695        cost_ref_traj = 1.8549259545339369e3
696        u_norm_ref_traj = 3.0455503580669716e2
697        cost_ref = 1.858428662785409e3
698        u_norm_ref = 3.0507636243132043e2
699       
700        # Get nominal and initial trajectories
701        opts = self.optimize_options(op, self.algorithm)
702        opts['n_e'] = 40
703        opts['n_cp'] = 2
704        res = op.optimize(self.algorithm, opts)
705        assert_results(res, cost_ref_traj, u_norm_ref_traj)
706        try:
707            os.remove("cstr_nom_traj_result.txt")
708        except OSError:
709            pass
710        os.rename("CSTR_CSTR_Opt_Bounds_Lagrange_result.txt",
711                  "cstr_nom_traj_result.txt")
712       
713        # Optimize using only initial trajectories
714        n_e = 75
715        n_cp = 4
716        opts['n_e'] = n_e
717        opts['n_cp'] = n_cp
718        opts['init_traj'] = ResultDymolaTextual("cstr_nom_traj_result.txt")
719        res = op.optimize(self.algorithm, opts)
720        assert_results(res, cost_ref, u_norm_ref)
721       
722        # Optimize using nominal and initial trajectories
723        opts['nominal_traj'] = ResultDymolaTextual("cstr_nom_traj_result.txt")
724        opts['nominal_traj_mode'] = {'_default_mode': 'time-variant'}
725        res = op.optimize(self.algorithm, opts)
726        assert_results(res, cost_ref, u_norm_ref)
727        col = res.solver
728        xx_init = col.get_xx_init()
729        N.testing.assert_equal(sum(xx_init == 1.),
730                               (n_e * n_cp + 1) * 4 + 2 * n_e)
731
732        # #3509
733        #~ # Test with eliminated continuity variables
734        #~ opts['eliminate_cont_var'] = True
735        #~ res = op.optimize(self.algorithm, opts)
736        #~ assert_results(res, cost_ref, u_norm_ref)
737        #~
738        #~ # Test with eliminated continuity and derivative variables
739        #~ opts['eliminate_der_var'] = True
740        #~ res = op.optimize(self.algorithm, opts)
741        #~ assert_results(res, cost_ref, u_norm_ref)
742       
743        # Test disabling scaling
744        opts['variable_scaling'] = False
745        res = op.optimize(self.algorithm, opts)
746        assert_results(res, cost_ref, u_norm_ref)
747       
748    @testattr(casadi_base = True)
749    def test_nominal_traj_cstr_lg(self):
750        """Test optimizing a CSTR using nominal and initial trajectories."""
751        op = self.cstr_lagrange_op
752       
753        # References values
754        #cost_ref_traj = 1.8549259545339369e3
755        cost_ref_traj = 1.86034214743738e3
756        #u_norm_ref_traj = 3.0455503580669716e2
757        u_norm_ref_traj = 3.053311192625306e2
758        cost_ref = 1.858428662785409e3
759        cost_ref_nominal = 1.6584308942977912e3
760        #u_norm_ref = 3.0507636243132043e2
761        u_norm_ref = 3.049332451630398e2
762        u_norm_ref_nominal = 3.0106051481333594e2
763       
764        # Get nominal and initial trajectories
765        opts = self.optimize_options(op, self.algorithm)
766        opts['n_e'] = 40
767        opts['n_cp'] = 2
768        opts["discr"] = "LG"
769        res = op.optimize(self.algorithm, opts)
770        assert_results(res, cost_ref_traj, u_norm_ref_traj)
771        try:
772            os.remove("cstr_nom_traj_result.txt")
773        except OSError:
774            pass
775        os.rename("CSTR_CSTR_Opt_Bounds_Lagrange_result.txt",
776                  "cstr_nom_traj_result.txt")
777       
778        # Optimize using only initial trajectories
779        n_e = 75
780        n_cp = 4
781        opts['n_e'] = n_e
782        opts['n_cp'] = n_cp
783        opts['init_traj'] = ResultDymolaTextual("cstr_nom_traj_result.txt")
784        res = op.optimize(self.algorithm, opts)
785        assert_results(res, cost_ref, u_norm_ref)
786       
787        # Optimize using nominal and initial trajectories
788        opts['nominal_traj'] = ResultDymolaTextual("cstr_nom_traj_result.txt")
789        opts['nominal_traj_mode'] = {'_default_mode': 'time-variant'}
790        res = op.optimize(self.algorithm, opts)
791        assert_results(res, cost_ref_nominal, u_norm_ref_nominal)
792        col = res.solver
793        xx_init = col.get_xx_init()
794        #N.testing.assert_equal(sum(xx_init == 1.),
795        #                       (n_e * n_cp + 1) * 4 + 2 * n_e)
796
797        # #3509
798        #~ # Test with eliminated continuity variables
799        #~ opts['eliminate_cont_var'] = True
800        #~ res = op.optimize(self.algorithm, opts)
801        #~ assert_results(res, cost_ref, u_norm_ref)
802        #~
803        #~ # Test with eliminated continuity and derivative variables
804        #~ opts['eliminate_der_var'] = True
805        #~ res = op.optimize(self.algorithm, opts)
806        #~ assert_results(res, cost_ref, u_norm_ref)
807       
808        # Test disabling scaling
809        opts['variable_scaling'] = False
810        res = op.optimize(self.algorithm, opts)
811        assert_results(res, cost_ref, u_norm_ref)
812
813    @testattr(casadi_base = True)
814    def test_nominal_traj_mode(self):
815        """Test nominal_traj_mode on the CSTR."""
816        op = self.cstr_lagrange_op
817       
818        # References values
819        cost_ref = 1.8549259545339369e3
820        u_norm_ref = 3.0455503580669716e2
821       
822        # Get nominal and initial trajectories
823        opts = self.optimize_options(op, self.algorithm)
824        opts['n_e'] = 40
825        opts['n_cp'] = 2
826        res = op.optimize(self.algorithm, opts)
827        assert_results(res, cost_ref, u_norm_ref)
828        try:
829            os.remove("cstr_nom_traj_result.txt")
830        except OSError:
831            pass
832        os.rename("CSTR_CSTR_Opt_Bounds_Lagrange_result.txt",
833                  "cstr_nom_traj_result.txt")
834       
835        # Time-variant
836        opts['nominal_traj'] = ResultDymolaTextual("cstr_nom_traj_result.txt")
837        opts['nominal_traj_mode'] = {'_default_mode': 'time-variant'}
838        res = op.optimize(self.algorithm, opts)
839        assert_results(res, cost_ref, u_norm_ref)
840       
841        # Affine
842        opts['nominal_traj_mode'] = {'_default_mode': 'affine'}
843        res = op.optimize(self.algorithm, opts)
844        assert_results(res, cost_ref, u_norm_ref)
845       
846        # Linear
847        opts['nominal_traj_mode'] = {'_default_mode': 'linear'}
848        res = op.optimize(self.algorithm, opts)
849        assert_results(res, cost_ref, u_norm_ref)
850       
851        # Attribute
852        opts['nominal_traj_mode'] = {'_default_mode': 'attribute'}
853        res = op.optimize(self.algorithm, opts)
854        assert_results(res, cost_ref, u_norm_ref)
855       
856        # Check impossible time-variant scaling
857        opts['nominal_traj_mode'] = {'_default_mode': 'affine',
858                                     'der(cstr.c)': 'time-variant'}
859        N.testing.assert_raises(CasadiCollocatorException, op.optimize,
860                                self.algorithm, opts)
861       
862        # Check changing one variable
863        opts['nominal_traj_mode'] = {'_default_mode': 'affine',
864                                     'cstr.c': 'time-variant'}
865        res = op.optimize(self.algorithm, opts)
866        assert_results(res, cost_ref, u_norm_ref)
867       
868        # Check alias variables
869        opts['nominal_traj_mode'] = {'_default_mode': 'affine',
870                                     'cstr.Tc': 'time-variant'}
871        res = op.optimize(self.algorithm, opts)
872        assert_results(res, cost_ref, u_norm_ref)
873        opts['nominal_traj_mode'] = {'_default_mode': 'affine',
874                                     'u': 'time-variant'}
875        res = op.optimize(self.algorithm, opts)
876        assert_results(res, cost_ref, u_norm_ref)
877        opts['nominal_traj_mode'] = {'_default_mode': 'affine',
878                                     'asdf': 'time-variant'}
879        N.testing.assert_raises(ValueError, op.optimize, self.algorithm,
880                                opts)
881
882    @testattr(casadi_base = True)
883    def test_cstr(self):
884        """
885        Test optimizing the CSTR.
886       
887        Tests both a Mayer cost with Gauss collocation and a Lagrange cost with
888        Radau collocation.
889        """
890        mayer_op = self.cstr_mayer_op
891        lagrange_op = self.cstr_lagrange_op
892       
893        # References values
894        cost_ref = 1.8576873858261e3
895        u_norm_ref = 3.050971000653911e2
896       
897        # Mayer
898        opts = self.optimize_options(mayer_op, self.algorithm)
899        opts['discr'] = "LG"
900        res = mayer_op.optimize(self.algorithm, opts)
901        assert_results(res, cost_ref, u_norm_ref)
902       
903        # Lagrange
904        opts['discr'] = "LGR"
905        res = lagrange_op.optimize(self.algorithm, opts)
906        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=5e-3)
907       
908    @testattr(casadi_base = True)
909    def test_cstr_checkpoint(self):
910        """
911        Test optimizing the CSTR with check point.
912       
913        Tests both a Mayer cost with Gauss collocation and a Lagrange cost with
914        Radau collocation.
915        """
916        mayer_op = self.cstr_mayer_op
917        lagrange_op = self.cstr_lagrange_op
918       
919        # References values
920        cost_ref = 1.8576873858261e3
921        u_norm_ref = 3.050971000653911e2
922       
923        # Mayer
924        opts = self.optimize_options(mayer_op, self.algorithm)
925        opts['discr'] = "LG"
926        opts['checkpoint'] = True
927        res = mayer_op.optimize(self.algorithm, opts)
928        assert_results(res, cost_ref, u_norm_ref)
929       
930        # Lagrange
931        opts['discr'] = "LGR"
932        opts['checkpoint'] = True
933        res = lagrange_op.optimize(self.algorithm, opts)
934        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=5e-3)   
935
936    @testattr(casadi_base = True)
937    def test_parameter_estimation(self):
938        """
939        Test a parameter estimation example with and without scaling.
940       
941        WARNING: This test is very slow when using IPOPT with the linear solver
942        MUMPS.
943        """
944        op = self.second_order_par_est_op
945       
946        # Reference values
947        w_ref = 1.048589
948        z_ref = 0.470934
949       
950        # Measurements
951        y_meas = N.array([0.01463904, 0.35424225, 0.94776816, 1.20116167,
952                          1.17283905, 1.03631145, 1.0549561, 0.94827652,
953                          1.0317119, 1.04010453, 1.08012155])
954        t_meas = N.linspace(0., 10., num=len(y_meas))
955        data = N.vstack([t_meas, y_meas])
956       
957        # Measurement data
958        Q = N.array([[1.]])
959        quad_pen = OrderedDict()
960        quad_pen['y'] = data
961        external_data = ExternalData(quad_pen=quad_pen, Q=Q)
962       
963        # Optimize without scaling
964        opts = self.optimize_options(op, self.algorithm)
965        opts['external_data'] = external_data
966        opts['variable_scaling'] = False
967        opts['n_e'] = 16
968        res = op.optimize(self.algorithm, opts)
969       
970        w_unscaled = res.final('w')
971        z_unscaled = res.final('z')
972        N.testing.assert_allclose(w_unscaled, w_ref, 1e-2)
973        N.testing.assert_allclose(z_unscaled, z_ref, 1e-2)
974       
975        # Optimize with scaling
976        opts['variable_scaling'] = True
977        res = op.optimize(self.algorithm, opts)
978        w_scaled = res.final('w')
979        z_scaled = res.final('z')
980        N.testing.assert_allclose(w_scaled, w_ref, 1e-2)
981        N.testing.assert_allclose(z_scaled, z_ref, 1e-2)
982
983    @testattr(casadi_base = True)
984    def test_parameter_estimation_traj(self):
985        """
986        Test estimation with and without initial and nominal trajectories.
987        """
988        model = self.second_order_model
989        model.reset()
990        op = self.second_order_par_est_op
991       
992        # Simulate with initial guesses
993        model.set('w', 1.3)
994        model.set('z', 0.3)
995        sim_res = model.simulate(final_time=15.)
996       
997        # Reference values
998        w_ref = 1.048589
999        z_ref = 0.470934
1000       
1001        # Measurements
1002        y_meas = N.array([0.01463904, 0.35424225, 0.94776816, 1.20116167,
1003                          1.17283905, 1.03631145, 1.0549561, 0.94827652,
1004                          1.0317119, 1.04010453, 1.08012155])
1005        t_meas = N.linspace(0., 10., num=len(y_meas))
1006        data = N.vstack([t_meas, y_meas])
1007       
1008        # External data
1009        Q = N.array([[1.]])
1010        quad_pen = OrderedDict()
1011        quad_pen['y'] = data
1012        external_data = ExternalData(quad_pen=quad_pen, Q=Q)
1013       
1014        # Optimize without scaling
1015        opts = self.optimize_options(op, self.algorithm)
1016        opts['external_data'] = external_data
1017        opts['n_e'] = 16
1018        opt_res = op.optimize(self.algorithm, opts)
1019       
1020        # Assert results
1021        w = opt_res.final('w')
1022        z = opt_res.final('z')
1023        N.testing.assert_allclose(w, w_ref, 1e-2)
1024        N.testing.assert_allclose(z, z_ref, 1e-2)
1025       
1026        # Optimize with trajectories
1027        opts['init_traj'] = sim_res
1028        opts['nominal_traj'] = sim_res
1029        traj_res = op.optimize(self.algorithm, opts)
1030        w_traj = traj_res.final('w')
1031        z_traj = traj_res.final('z')
1032        N.testing.assert_allclose(w_traj, w_ref, 1e-2)
1033        N.testing.assert_allclose(z_traj, z_ref, 1e-2)
1034
1035    @testattr(casadi_base = True)
1036    def test_qt_par_est_quad_pen(self):
1037        """
1038        Test parameter estimation for the quad tank with quad_pen inputs.
1039        """
1040        data = loadmat(path_to_data + '/qt_par_est_data.mat', appendmat=False)
1041        model = self.qt_model
1042        model.reset()
1043        op = self.qt_par_est_op
1044        a_ref = [0.02656702, 0.02713898]
1045       
1046        # Extract data series
1047        t_meas = data['t'][6000::100, 0] - 60
1048        y1_meas = data['y1_f'][6000::100, 0]/100
1049        y2_meas = data['y2_f'][6000::100, 0]/100
1050        y3_meas = data['y3_d'][6000::100, 0]/100
1051        y4_meas = data['y4_d'][6000::100, 0]/100
1052        u1 = data['u1_d'][6000::100, 0]
1053        u2 = data['u2_d'][6000::100, 0]
1054       
1055        # Simulate
1056        u = N.transpose(N.vstack((t_meas, u1, u2)))
1057        sim_res = model.simulate(input=(['u1', 'u2'], u), start_time=0.,
1058                                 final_time=60.)
1059       
1060        # Create external data
1061        Q = N.diag([1., 1., 10., 10.])
1062        data_x1 = N.vstack([t_meas, y1_meas])
1063        data_x2 = N.vstack([t_meas, y2_meas])
1064        data_u1 = N.vstack([t_meas, u1])
1065        data_u2 = N.vstack([t_meas, u2])
1066        quad_pen = OrderedDict()
1067        quad_pen['qt.x1'] = data_x1
1068        quad_pen['qt.x2'] = data_x2
1069        quad_pen['u1'] = data_u1
1070        quad_pen['u2'] = data_u2
1071        external_data = ExternalData(Q=Q, quad_pen=quad_pen)
1072       
1073        # Unconstrained
1074        opts = self.optimize_options(op, self.algorithm)
1075        opts['n_e'] = 60
1076        opts['init_traj'] = sim_res
1077        opts['external_data'] = external_data
1078        opt_res = op.optimize(self.algorithm, options=opts)
1079        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1080                                                 opt_res.final("qt.a2")]),
1081                                  a_ref, rtol=1e-4)
1082       
1083        # Unconstrained with nominal trajectories
1084        opts['nominal_traj'] = sim_res
1085        opt_res = op.optimize(self.algorithm, options=opts)
1086        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1087                                                 opt_res.final("qt.a2")]),
1088                                  a_ref, rtol=1e-4)
1089
1090    @testattr(casadi_base = True)
1091    def test_qt_par_est_eliminated(self):
1092        """
1093        Test parameter estimation for the quad tank with eliminated inputs.
1094        """
1095        data = loadmat(path_to_data + '/qt_par_est_data.mat', appendmat=False)
1096        model = self.qt_model
1097        model.reset()
1098        op = self.qt_par_est_op
1099        op_2 = self.qt_par_est_degenerate_op
1100        a_ref = [0.02656702, 0.02713898]
1101       
1102        # Extract data series
1103        t_meas = data['t'][6000::100, 0] - 60
1104        y1_meas = data['y1_f'][6000::100, 0]/100
1105        y2_meas = data['y2_f'][6000::100, 0]/100
1106        y3_meas = data['y3_d'][6000::100, 0]/100
1107        y4_meas = data['y4_d'][6000::100, 0]/100
1108        u1 = data['u1_d'][6000::100, 0]
1109        u2 = data['u2_d'][6000::100, 0]
1110       
1111        # Simulate
1112        u = N.transpose(N.vstack((t_meas, u1, u2)))
1113        sim_res = model.simulate(input=(['u1', 'u2'], u), start_time=0.,
1114                                 final_time=60.)
1115       
1116        # Create external data
1117        Q = N.diag([1., 1.])
1118        data_x1 = N.vstack([t_meas, y1_meas])
1119        data_x2 = N.vstack([t_meas, y2_meas])
1120        data_u1 = N.vstack([t_meas, u1])
1121        data_u2 = N.vstack([t_meas, u2])
1122        quad_pen = OrderedDict()
1123        quad_pen['qt.x1'] = data_x1
1124        quad_pen['qt.x2'] = data_x2
1125        eliminated = OrderedDict()
1126        eliminated['u1'] = data_u1
1127        eliminated['u2'] = data_u2
1128        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1129                                     eliminated=eliminated)
1130       
1131        # Eliminated
1132        opts = self.optimize_options(op, self.algorithm)
1133        opts['n_e'] = 30
1134        opts['init_traj'] = sim_res
1135        opts['external_data'] = external_data
1136        opt_res = op.optimize(self.algorithm, opts)
1137        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1138                                                 opt_res.final("qt.a2")]),
1139                                  a_ref, rtol=1e-3)
1140       
1141        # Eliminated with nominal trajectories
1142        opts['nominal_traj'] = sim_res
1143        opt_res = op.optimize(self.algorithm, opts)
1144        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1145                                                 opt_res.final("qt.a2")]),
1146                                  a_ref, rtol=1e-3)
1147       
1148        # Inconsistent bound on eliminated input
1149        op.getVariable('u1').setMin(5.1)
1150        N.testing.assert_raises(CasadiCollocatorException,
1151                                op.optimize, self.algorithm, opts)
1152        op.getVariable('u1').setMin(-N.inf)
1153       
1154        # Point constraint on eliminated input
1155        opts2 = self.optimize_options(op_2, self.algorithm)
1156        opts2['n_e'] = 30
1157        opts2['init_traj'] = sim_res
1158        opts2['external_data'] = external_data
1159        N.testing.assert_raises(CasadiCollocatorException,
1160                                op_2.optimize, self.algorithm, opts2)
1161       
1162        # Eliminate state
1163        Q = N.diag([1.])
1164        quad_pen = OrderedDict()
1165        quad_pen['qt.x1'] = data_x1
1166        eliminated = OrderedDict()
1167        eliminated['u1'] = data_u1
1168        eliminated['u2'] = data_u2
1169        eliminated['qt.x2'] = data_x2
1170        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1171                                     eliminated=eliminated)
1172        opts['external_data'] = external_data
1173        N.testing.assert_raises(CasadiCollocatorException,
1174                                op.optimize, self.algorithm, opts)
1175
1176        # Eliminate non-existing variable
1177        del eliminated['qt.x2']
1178        eliminated['does_not_exist'] = data_x2
1179        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1180                                     eliminated=eliminated)
1181        opts['external_data'] = external_data
1182        N.testing.assert_raises(CasadiCollocatorException,
1183                                op.optimize, self.algorithm, opts)
1184
1185    @testattr(casadi_base = True)
1186    def test_qt_par_est_constr_quad_pen(self):
1187        """
1188        Test parameter estimation for quad tank with constr_quad_pen inputs.
1189        """
1190        data = loadmat(path_to_data + '/qt_par_est_data.mat', appendmat=False)
1191        model = self.qt_model
1192        model.reset()
1193        op = self.qt_par_est_op
1194        a_ref = [0.02656702, 0.02713898]
1195       
1196        # Extract data series
1197        t_meas = data['t'][6000::100, 0] - 60
1198        y1_meas = data['y1_f'][6000::100, 0]/100
1199        y2_meas = data['y2_f'][6000::100, 0]/100
1200        y3_meas = data['y3_d'][6000::100, 0]/100
1201        y4_meas = data['y4_d'][6000::100, 0]/100
1202        u1 = data['u1_d'][6000::100, 0]
1203        u2 = data['u2_d'][6000::100, 0]
1204       
1205        # Simulate
1206        u = N.transpose(N.vstack((t_meas, u1, u2)))
1207        sim_res = model.simulate(input=(['u1', 'u2'], u), start_time=0.,
1208                                 final_time=60.)
1209       
1210        # Create external data
1211        Q = N.diag([1., 1., 10., 10.])
1212        data_x1 = N.vstack([t_meas, y1_meas])
1213        data_x2 = N.vstack([t_meas, y2_meas])
1214        data_u1 = N.vstack([t_meas, u1])
1215        data_u2 = N.vstack([t_meas, u2])
1216        quad_pen = OrderedDict()
1217        quad_pen['qt.x1'] = data_x1
1218        quad_pen['qt.x2'] = data_x2
1219        constr_quad_pen = OrderedDict()
1220        constr_quad_pen['u1'] = data_u1
1221        constr_quad_pen['u2'] = data_u2
1222        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1223                                     constr_quad_pen=constr_quad_pen)
1224       
1225        # constr_quad_pen
1226        opts = self.optimize_options(op, self.algorithm)
1227        opts['n_e'] = 30
1228        opts['init_traj'] = sim_res
1229        opts['external_data'] = external_data
1230        opt_res = op.optimize(self.algorithm, opts)
1231        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1232                                                 opt_res.final("qt.a2")]),
1233                                  a_ref, rtol=1e-3)
1234       
1235        # constr_quad_pen with nominal trajectories
1236        opts['nominal_traj'] = sim_res
1237        opt_res = op.optimize(self.algorithm, opts)
1238        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1239                                                 opt_res.final("qt.a2")]),
1240                                  a_ref, rtol=1e-3)
1241       
1242        # Inconsistent bound on constrained input
1243        op.getVariable('u1').setMin(5.1)
1244        N.testing.assert_raises(CasadiCollocatorException,
1245                                op.optimize, self.algorithm, opts)
1246        op.getVariable('u1').setMin(-N.inf)
1247       
1248        # Constrain state
1249        quad_pen = OrderedDict()
1250        quad_pen['qt.x1'] = data_x1
1251        constr_quad_pen = OrderedDict()
1252        constr_quad_pen['u1'] = data_u1
1253        constr_quad_pen['u2'] = data_u2
1254        constr_quad_pen['qt.x2'] = data_x2
1255        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1256                                     constr_quad_pen=constr_quad_pen)
1257        opts['external_data'] = external_data
1258        N.testing.assert_raises(VariableNotFoundError,
1259                                op.optimize, self.algorithm, opts)
1260
1261    @testattr(casadi_base = True)
1262    def test_qt_par_est_semi_eliminated(self):
1263        """
1264        Test parameter estimation for the quad tank with 1 eliminated input.
1265        """
1266        data = loadmat(path_to_data + '/qt_par_est_data.mat', appendmat=False)
1267        model = self.qt_model
1268        model.reset()
1269        op = self.qt_par_est_op
1270        a_ref = [0.02656702, 0.02713898]
1271       
1272        # Extract data series
1273        t_meas = data['t'][6000::100, 0] - 60
1274        y1_meas = data['y1_f'][6000::100, 0]/100
1275        y2_meas = data['y2_f'][6000::100, 0]/100
1276        y3_meas = data['y3_d'][6000::100, 0]/100
1277        y4_meas = data['y4_d'][6000::100, 0]/100
1278        u1 = data['u1_d'][6000::100, 0]
1279        u2 = data['u2_d'][6000::100, 0]
1280       
1281        # Simulate
1282        u = N.transpose(N.vstack((t_meas, u1, u2)))
1283        sim_res = model.simulate(input=(['u1', 'u2'], u), start_time=0.,
1284                                     final_time=60.)
1285       
1286        # Create external data
1287        Q = N.diag([1., 1., 10.])
1288        data_x1 = N.vstack([t_meas, y1_meas])
1289        data_x2 = N.vstack([t_meas, y2_meas])
1290        data_u1 = N.vstack([t_meas, u1])
1291        data_u2 = N.vstack([t_meas, u2])
1292        quad_pen = OrderedDict()
1293        quad_pen['qt.x1'] = data_x1
1294        quad_pen['qt.x2'] = data_x2
1295        quad_pen['u1'] = data_u1
1296        eliminated = OrderedDict()
1297        eliminated['u2'] = data_u2
1298        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1299                                     eliminated=eliminated)
1300       
1301        # Eliminate u2
1302        opts = self.optimize_options(op, self.algorithm)
1303        opts['n_e'] = 60
1304        opts['init_traj'] = sim_res
1305        opts['external_data'] = external_data
1306        opt_res = op.optimize(self.algorithm, opts)
1307        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1308                                                 opt_res.final("qt.a2")]),
1309                                  a_ref, rtol=1e-4)
1310       
1311        # Eliminate u2 with nominal trajectories
1312        opts['nominal_traj'] = sim_res
1313        opt_res = op.optimize(self.algorithm, opts)
1314        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1315                                                 opt_res.final("qt.a2")]),
1316                                  a_ref, rtol=1e-4)
1317       
1318        # Eliminate u1
1319        quad_pen = OrderedDict()
1320        quad_pen['qt.x1'] = data_x1
1321        quad_pen['qt.x2'] = data_x2
1322        quad_pen['u2'] = data_u2
1323        eliminated = OrderedDict()
1324        eliminated['u1'] = data_u1
1325        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1326                                     eliminated=eliminated)
1327        opts['nominal_traj'] = None
1328        opts['external_data'] = external_data
1329        opt_res = op.optimize(self.algorithm, opts)
1330        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1331                                                 opt_res.final("qt.a2")]),
1332                                  a_ref, rtol=1e-4)
1333       
1334        # Eliminate u1 with nominal trajectories
1335        opts['nominal_traj'] = sim_res
1336        opt_res = op.optimize(self.algorithm, opts)
1337        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1338                                                 opt_res.final("qt.a2")]),
1339                                  a_ref, rtol=1e-4)
1340
1341    @testattr(casadi_base = True)
1342    def test_qt_par_est_user_interpolator(self):
1343        """
1344        Test parameter estimation for the quad tank with user-defined function.
1345        """
1346        data = loadmat(path_to_data + '/qt_par_est_data.mat', appendmat=False)
1347        model = self.qt_model
1348        model.reset()
1349        op = self.qt_par_est_op
1350        a_ref = [0.02656702, 0.02713898]
1351       
1352        # Extract data series
1353        t_meas = data['t'][6000::100, 0] - 60
1354        y1_meas = data['y1_f'][6000::100, 0]/100
1355        y2_meas = data['y2_f'][6000::100, 0]/100
1356        y3_meas = data['y3_d'][6000::100, 0]/100
1357        y4_meas = data['y4_d'][6000::100, 0]/100
1358        u1 = data['u1_d'][6000::100, 0]
1359        u2 = data['u2_d'][6000::100, 0]
1360       
1361        # Simulate
1362        u = N.transpose(N.vstack((t_meas, u1, u2)))
1363        sim_res = model.simulate(input=(['u1', 'u2'], u), start_time=0.,
1364                                     final_time=60.)
1365       
1366        # Create external data
1367        Q = N.diag([1., 1., 10.])
1368        data_x1 = N.vstack([t_meas, y1_meas])
1369        data_x2 = N.vstack([t_meas, y2_meas])
1370        data_u2 = N.vstack([t_meas, u2])
1371        linear_u1 = TrajectoryLinearInterpolation(t_meas,
1372                                                  u1.reshape([-1, 1]))
1373        user_u1 = linear_u1.eval
1374        quad_pen = OrderedDict()
1375        quad_pen['qt.x1'] = data_x1
1376        quad_pen['qt.x2'] = data_x2
1377        quad_pen['u2'] = data_u2
1378        eliminated = OrderedDict()
1379        eliminated['u1'] = user_u1
1380        external_data = ExternalData(Q=Q, quad_pen=quad_pen,
1381                                     eliminated=eliminated)
1382       
1383        # Semi-eliminated with user-defined interpolation function
1384        opts = self.optimize_options(op, self.algorithm)
1385        opts['n_e'] = 60
1386        opts['external_data'] = external_data
1387        opt_res = op.optimize(self.algorithm, opts)
1388        N.testing.assert_allclose(1e4 * N.array([opt_res.final("qt.a1"),
1389                                                 opt_res.final("qt.a2")]),
1390                                  a_ref, rtol=1e-4)
1391       
1392        # Inconsistent bounds on eliminated input with user-defined function
1393        op.getVariable('u1').setMin(5.1)
1394        N.testing.assert_raises(CasadiCollocatorException,
1395                                op.optimize, self.algorithm, opts)
1396        op.getVariable('u1').setMin(-N.inf)
1397
1398    @testattr(casadi_base = True)
1399    def test_vdp_minimum_time(self):
1400        """
1401        Test solving minimum time problems based on the VDP oscillator.
1402       
1403        Tests one problem where the time is manually scaled. Tests two where
1404        the time is automatically scaled by the compiler, where one of them
1405        starts at time 0 and the other starts at time 5.
1406        """
1407        scaled_op = self.vdp_scaled_min_time_op
1408        unscaled_op = self.vdp_unscaled_min_time_op
1409        nonzero_start_op = \
1410                self.vdp_min_time_nonzero_start_op
1411
1412        # References values
1413        cost_ref = 2.2811590707107996e0
1414        u_norm_ref = 9.991517452037317e-1
1415       
1416        # Scaled, Radau
1417        opts = self.optimize_options(scaled_op, self.algorithm)
1418        opts['discr'] = "LGR"
1419        res = scaled_op.optimize(self.algorithm, opts)
1420        assert_results(res, cost_ref, u_norm_ref)
1421       
1422        # Scaled, Gauss
1423        opts['discr'] = "LG"
1424        res = scaled_op.optimize(self.algorithm, opts)
1425        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1426       
1427        # Unscaled, Radau
1428        opts['discr'] = "LGR"
1429        res = unscaled_op.optimize(self.algorithm, opts)
1430        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1431        res_radau = res
1432       
1433        # Unscaled, Gauss
1434        opts['discr'] = "LG"
1435        res = unscaled_op.optimize(self.algorithm, opts)
1436        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1437
1438        # Unscaled, non-zero start, Radau
1439        opts['discr'] = "LGR"
1440        res = nonzero_start_op.optimize(self.algorithm, opts)
1441        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1442        N.testing.assert_allclose(res['time'][[0, -1]], [5., 7.28115907],
1443                                  rtol=5e-3)
1444       
1445        # Unscaled, non-zero start, Gauss
1446        opts['discr'] = "LG"
1447        res = nonzero_start_op.optimize(self.algorithm, opts)
1448        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1449        N.testing.assert_allclose(res['time'][[0, -1]], [5., 7.28128126],
1450                                  rtol=5e-3)
1451
1452        # Unscaled, nominal trajectories, Radau
1453        opts['discr'] = "LGR"
1454        opts['nominal_traj'] = res_radau
1455        res = unscaled_op.optimize(self.algorithm, opts)
1456        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2)
1457
1458    @testattr(casadi_base = True)
1459    def test_reorder(self):
1460        """Test reordered NLP."""
1461        op = self.cstr_lagrange_op
1462       
1463        # References values
1464        cost_ref = 1.858428662785409e3
1465        u_norm_ref = 305.5673000177923
1466       
1467        # Default order
1468        opts = self.optimize_options(op, self.algorithm)
1469        opts['write_scaled_result'] = True
1470        opts['order'] = "default"
1471        res = op.optimize(self.algorithm, opts)
1472        assert_results(res, cost_ref, u_norm_ref)
1473       
1474        # Reverse order
1475        opts['order'] = "reverse"
1476        res = op.optimize(self.algorithm, opts)
1477        assert_results(res, cost_ref, u_norm_ref)
1478       
1479        # Random order #1
1480        N.random.seed(1)
1481        opts['order'] = "random"
1482        res = op.optimize(self.algorithm, opts)
1483        assert_results(res, cost_ref, u_norm_ref)
1484       
1485        # Random order #2
1486        res = op.optimize(self.algorithm, opts)
1487        assert_results(res, cost_ref, u_norm_ref)
1488
1489    @testattr(casadi_base = True)
1490    def test_cstr_minimum_time(self):
1491        """
1492        Test solving minimum time problems based on the CSTR.
1493       
1494        Tests both a problem where the time is manually scaled, and one where
1495        the time is automatically scaled by the compiler.
1496        """
1497        scaled_op = self.cstr_scaled_min_time_op
1498        unscaled_op = self.cstr_unscaled_min_time_op
1499       
1500        # References values
1501        cost_ref = 1.1637020114180874e2
1502        u_norm_ref = 3.0668821961641106e2
1503       
1504        # Scaled, Radau
1505        opts = self.optimize_options(scaled_op, self.algorithm)
1506        opts['discr'] = "LGR"
1507        opts['IPOPT_options']['linear_solver'] = "mumps"
1508        res = scaled_op.optimize(self.algorithm, opts)
1509        assert_results(res, cost_ref, u_norm_ref, input_name="Tc")
1510       
1511        # Scaled, Gauss
1512        opts['discr'] = "LG"
1513        opts['IPOPT_options']['linear_solver'] = "mumps"
1514        res = scaled_op.optimize(self.algorithm, opts)
1515        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2,
1516                       input_name="Tc")
1517       
1518        # Unscaled, Radau
1519        opts['discr'] = "LGR"
1520        opts['IPOPT_options']['linear_solver'] = "mumps"
1521        res = unscaled_op.optimize(self.algorithm, opts)
1522        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2,
1523                       input_name="Tc")
1524       
1525        # Unscaled, Gauss
1526        opts['discr'] = "LG"
1527        opts['IPOPT_options']['linear_solver'] = "mumps"
1528        res = unscaled_op.optimize(self.algorithm, opts)
1529        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=1e-2,
1530                       input_name="Tc")
1531
1532    @testattr(casadi_base = True)
1533    def test_path_constraints(self):
1534        """Test a simple path constraint with and without exact Hessian."""
1535        op = self.vdp_constraints_mayer_op
1536       
1537        # References values
1538        cost_ref = 5.273481330869811e0
1539        u_norm_ref = 3.2936323844551e-1
1540       
1541        # Without exact Hessian
1542        opts = self.optimize_options(op, self.algorithm)
1543        opts['IPOPT_options']['hessian_approximation'] = "limited-memory"
1544        res = op.optimize(self.algorithm, opts)
1545        assert_results(res, cost_ref, u_norm_ref)
1546       
1547        # With exact Hessian
1548        opts['IPOPT_options']['hessian_approximation'] = "exact"
1549        res = op.optimize(self.algorithm, opts)
1550        assert_results(res, cost_ref, u_norm_ref)
1551
1552    @testattr(casadi_base = True)
1553    def test_initial_equations(self):
1554        """Test initial equations with and without eliminated derivatives."""
1555        op = self.vdp_initial_equations_op
1556       
1557        # References values
1558        cost_ref = 4.7533158101416788e0
1559        u_norm_ref = 5.18716394291585e-1
1560       
1561        # Without derivative elimination
1562        opts = self.optimize_options(op, self.algorithm)
1563        opts['eliminate_der_var'] = False
1564        res = op.optimize(self.algorithm, opts)
1565        assert_results(res, cost_ref, u_norm_ref)
1566
1567        # #3509
1568        #~ # With derivative elimination
1569        #~ opts['eliminate_der_var'] = True
1570        #~ res = op.optimize(self.algorithm, opts)
1571        #~ assert_results(res, cost_ref, u_norm_ref)
1572
1573    @testattr(casadi_base = True)
1574    def test_element_lengths(self):
1575        """Test non-uniformly distributed elements."""
1576        op = self.vdp_bounds_mayer_op
1577       
1578        opts = self.optimize_options(op, self.algorithm)
1579        opts['n_e'] = 23
1580        opts['hs'] = N.array(4 * [0.01] + 2 * [0.05] + 10 * [0.02] +
1581                             5 * [0.02] + 2 * [0.28])
1582        res = op.optimize(self.algorithm, opts)
1583        assert_results(res, 3.174936706809e0, 3.707273799325e-1)
1584
1585    @testattr(casadi_base = True)
1586    def test_free_element_lengths(self):
1587        """Test optimized element lengths with both result modes."""
1588        op = self.vdp_bounds_lagrange_op
1589       
1590        # References values
1591        cost_ref = 3.3821187315826737e0
1592        u_norm_ref = 4.011979950965081e-1
1593       
1594        # Free element lengths data
1595        c = 0.5
1596        Q = N.eye(2)
1597        bounds = (0.5, 2.0)
1598        free_ele_data = FreeElementLengthsData(c, Q, bounds)
1599       
1600        # Set options shared by both result modes
1601        opts = self.optimize_options(op, self.algorithm)
1602        opts['n_e'] = 20
1603        opts['hs'] = "free"
1604        opts['free_element_lengths_data'] = free_ele_data
1605       
1606        # Collocation points
1607        opts['result_mode'] = "collocation_points"
1608        res = op.optimize(self.algorithm, opts)
1609        assert_results(res, cost_ref, u_norm_ref)
1610        indices = list(range(1, 4)) + list(range(opts['n_e'] - 3, opts['n_e']))
1611        values = N.array([0.5, 0.5, 0.5, 2.0, 2.0, 2.0])
1612        N.testing.assert_allclose(20. * res.h_opt[indices], values, 5e-3)
1613       
1614        # Element interpolation
1615        opts['result_mode'] = "element_interpolation"
1616        res = op.optimize(self.algorithm, opts)
1617        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=3e-2)
1618
1619    @testattr(casadi_base = True)
1620    def test_named_vars(self):
1621        """
1622        Test variable renaming.
1623
1624        This test is by no means thorough.
1625        """
1626        op = self.vdp_bounds_mayer_op
1627       
1628        # References values
1629        cost_ref = 1.353983656973385e0
1630        u_norm_ref = 2.4636859805244668e-1
1631
1632        # Common options
1633        opts = self.optimize_options(op, self.algorithm)
1634        opts['n_e'] = 2
1635
1636        # Without naming
1637        opts['named_vars'] = False
1638        res = op.optimize(self.algorithm, opts)
1639        assert_results(res, cost_ref, u_norm_ref)
1640
1641        # With renaming
1642        opts['named_vars'] = True
1643        res_renaming = op.optimize(self.algorithm, opts)
1644        assert_results(res_renaming, cost_ref, u_norm_ref)
1645
1646        # Compare tenth and twentieth equality constraints
1647        c_e_10 = res_renaming.solver.get_equality_constraint()[self.test_named_vars_c_e_10_id()]
1648        c_e_20 = res_renaming.solver.get_equality_constraint()[20]
1649       
1650        #comp_str1 = "SX((((der(x1)_1_1_d_sf*der(x1)_1_1)+der(x1)_1_1_e_sf)-((((1-sq(((x2_1_1_d_sf*x2_1_1)+x2_1_1_e_sf)))*((der(x2)_1_1_d_sf*der(x2)_1_1)+der(x2)_1_1_e_sf))-((x2_1_1_d_sf*x2_1_1)+x2_1_1_e_sf))+((u_1_1_d_sf*u_1_1)+u_1_1_e_sf))))"
1651        #comp_str1 = "SX((((der(x1)_d_sf*der(x1)_1_1)+der(x1)_e_sf)-((((1-sq(((x2_d_sf*x2_1_1)+x2_e_sf)))*((der(x2)_d_sf*der(x2)_1_1)+der(x2)_e_sf))-((x2_d_sf*x2_1_1)+x2_e_sf))+((u_d_sf*u_1_1)+u_e_sf))))"
1652        comp_str1 = "SX((der(x1)_1_1-((((1-sq(x2_1_1))*der(x2)_1_1)-x2_1_1)+u_1_1)))"
1653       
1654        #comp_str2 = "SX((((((-3*((x2_1_0_d_sf*x2_1_0)+x2_1_0_e_sf))+(5.53197*((x2_1_1_d_sf*x2_1_1)+x2_1_1_e_sf)))+(-7.53197*((x2_1_2_d_sf*x2_1_2)+x2_1_2_e_sf)))+(5*((x2_1_3_d_sf*x2_1_3)+x2_1_3_e_sf)))-(10*((der(x2)_1_3_d_sf*der(x2)_1_3)+der(x2)_1_3_e_sf))))"
1655        #comp_str2 = "SX((((((-3*((x2_d_sf*x2_1_0)+x2_e_sf))+(5.53197*((x2_d_sf*x2_1_1)+x2_e_sf)))+(-7.53197*((x2_d_sf*x2_1_2)+x2_e_sf)))+(5*((x2_d_sf*x2_1_3)+x2_e_sf)))-(10*((der(x2)_d_sf*der(x2)_1_3)+der(x2)_e_sf))))"
1656        comp_str2 = "SX((((((-3*x2_1_0)+(5.53197*x2_1_1))+(-7.53197*x2_1_2))+(5*x2_1_3))-(10*der(x2)_1_3)))"
1657       
1658        N.testing.assert_string_equal(
1659            repr(res_renaming.solver.get_named_var_expr(c_e_10)),
1660                 comp_str1)
1661        N.testing.assert_string_equal(
1662            repr(res_renaming.solver.get_named_var_expr(c_e_20)),
1663            comp_str2)
1664
1665    def test_named_vars_c_e_10_id(self):
1666        return 10;
1667   
1668    @testattr(casadi_base = True)
1669    def test_scaling(self):
1670        """
1671        Test optimizing the CSTR with and without scaling.
1672
1673        This test also tests writing both the unscaled and scaled result,
1674        eliminating derivative variables and setting nominal values
1675        post-compilation.
1676        """
1677        op = self.cstr_lagrange_op
1678       
1679        # References values
1680        cost_ref = 1.8576873858261e3
1681        u_norm_ref = 3.0556730059e2
1682       
1683        # Unscaled variables, with derivatives
1684        opts = self.optimize_options(op, self.algorithm)
1685        opts['variable_scaling'] = False
1686        opts['write_scaled_result'] = False
1687        opts['eliminate_der_var'] = False
1688        res = op.optimize(self.algorithm, opts)
1689        assert_results(res, cost_ref, u_norm_ref)
1690
1691        # 3509
1692        #~ # Scaled variables, unscaled result
1693        #~ # Eliminated derivatives
1694        #~ opts['variable_scaling'] = True
1695        #~ opts['write_scaled_result'] = False
1696        #~ opts['eliminate_der_var'] = True
1697        #~ res = op.optimize(self.algorithm, opts)
1698        #~ assert_results(res, cost_ref, u_norm_ref)
1699        #~ c_unscaled = res['cstr.c']
1700#~
1701        # Scaled variables, scaled result
1702        #~ # Eliminated derivatives
1703        #~ opts['variable_scaling'] = True
1704        #~ opts['write_scaled_result'] = True
1705        #~ opts['eliminate_der_var'] = True
1706        #~ res = op.optimize(self.algorithm, opts)
1707        #~ assert_results(res, cost_ref, u_norm_ref)
1708        #~ c_scaled = res['cstr.c']
1709        #~ N.testing.assert_allclose(c_unscaled, 1000. * c_scaled,
1710                                  #~ rtol=0, atol=1e-5)
1711        #~
1712        #~ # Scaled variables, scaled result, with updated nominal value
1713        #~ # Eliminated derivatives
1714        #~ op.set_nominal('cstr.c', 500.)
1715        #~ res = op.optimize(self.algorithm, opts)
1716        #~ op.set_nominal('cstr.c', 1000.)
1717        #~ assert_results(res, cost_ref, u_norm_ref)
1718        #~ c_scaled = res['cstr.c']
1719        #~ N.testing.assert_allclose(c_unscaled, 500. * c_scaled,
1720                                  #~ rtol=0, atol=1e-5)
1721
1722    @testattr(casadi_base = True)
1723    def test_result_file_name(self):
1724        """
1725        Test different result file names.
1726        """
1727        op = self.vdp_bounds_lagrange_op
1728       
1729        # Default file name
1730        try:
1731            os.remove("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt")
1732        except OSError:
1733            pass
1734        op.optimize(self.algorithm)
1735        assert(os.path.exists("VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt"))
1736       
1737        # Custom file name
1738        opts = self.optimize_options(op, self.algorithm)
1739        opts['result_file_name'] = "vdp_custom_file_name.txt"
1740        try:
1741            os.remove("vdp_custom_file_name.txt")
1742        except OSError:
1743            pass
1744        op.optimize(self.algorithm, opts)
1745        assert(os.path.exists("vdp_custom_file_name.txt"))
1746
1747    @testattr(casadi_base = True)
1748    def test_result_mode(self):
1749        """
1750        Test the two different result modes.
1751       
1752        The difference between the trajectories of the three result modes
1753        should be very small if n_e * n_cp is sufficiently large. Eliminating
1754        derivative variables is also tested.
1755        """
1756        op = self.vdp_bounds_lagrange_op
1757       
1758        # References values
1759        cost_ref = 3.17495094634053e0
1760        u_norm_ref = 2.84538299160e-1
1761       
1762        # Collocation points
1763        opts = self.optimize_options(op, self.algorithm)
1764        opts['n_e'] = 100
1765        opts['n_cp'] = 5
1766        opts['result_mode'] = "collocation_points"
1767        res = op.optimize(self.algorithm, opts)
1768        assert_results(res, cost_ref, u_norm_ref)
1769       
1770        # Element interpolation
1771        opts['result_mode'] = "element_interpolation"
1772        #~ opts['eliminate_der_var'] = True #3509
1773        opts['n_eval_points'] = 15
1774        res = op.optimize(self.algorithm, opts)
1775        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=5e-3)
1776       
1777        # Mesh points
1778        opts['result_mode'] = "mesh_points"
1779        opts['n_eval_points'] = 20 # Reset to default
1780        res = op.optimize(self.algorithm, opts)
1781        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=5e-3)
1782
1783    @testattr(casadi_base = True)
1784    def test_parameter_setting(self):
1785        """
1786        Test setting parameters post-compilation.
1787        """
1788        op = self.cstr_dependent_parameter_op
1789       
1790        # Reference values
1791        cost_ref_low = 1.2391821615924346e3
1792        u_norm_ref_low = 2.833724580055005e2
1793        cost_ref_default = 1.8576873858261e3
1794        u_norm_ref_default = 3.0556730059139556e2
1795       
1796        # Test lower EdivR
1797        op.set('cstr.EdivR', 8200)
1798        res_low = op.optimize(self.algorithm)
1799        assert_results(res_low, cost_ref_low, u_norm_ref_low)
1800       
1801        # Test default EdviR
1802        op.set('cstr.EdivR', 8750)
1803        res_default = op.optimize(self.algorithm)
1804        assert_results(res_default, cost_ref_default, u_norm_ref_default)
1805
1806        # Test dependent parameter setting
1807        N.testing.assert_raises(RuntimeError, op.set, 'cstr.F', 500)
1808
1809    @testattr(casadi_base = True)
1810    def test_blocking_factors(self):
1811        """Test blocking factors."""
1812        op = self.vdp_bounds_lagrange_op
1813
1814        # Check constant blocking factors
1815        opts = self.optimize_options(op, self.algorithm)
1816        opts['n_e'] = 40
1817        opts['n_cp'] = 3
1818        opts['blocking_factors'] = N.array(opts['n_e'] * [1])
1819        res = op.optimize(self.algorithm, opts)
1820        assert_results(res, 3.310907e0, 2.8718067e-1,
1821                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1822
1823        # Check varying blocking factors
1824        opts['n_e'] = 20
1825        opts['n_cp'] = 4
1826        opts['blocking_factors'] = [1, 1, 1, 1, 1, 2, 2, 2, 9]
1827        res = op.optimize(self.algorithm, opts)
1828        assert_results(res, 3.620908e0, 3.048898e-1,
1829                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1830
1831        # Check blocking factors with bound and penalty
1832        factors = {'u': [1, 1, 1, 1, 1, 2, 2, 2, 9]}
1833        du_quad_pen = {'u': 1}
1834        du_bounds = {'u': 0.6}
1835        bf = BlockingFactors(factors, du_quad_pen, du_bounds)
1836        opts['blocking_factors'] = bf
1837        res = op.optimize(self.algorithm, opts)
1838        assert_results(res, 3.883757e0, 3.068191e-1,
1839                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1840
1841    @testattr(casadi_base = True)
1842    def test_blocking_factors_cstr(self):
1843        """Test blocking factors for CSTR."""
1844        op = self.cstr_lagrange_op
1845        cost_ref = 1.873e3
1846        u_norm_ref = 3.053e2
1847        blocking_factors = [12, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1,
1848                            2, 2, 2, 3, 3, 6]
1849
1850        # Check blocking factors without scaling
1851        opts = self.optimize_options(op, self.algorithm)
1852        opts['blocking_factors'] = blocking_factors
1853        res = op.optimize(self.algorithm, opts)
1854        assert_results(res, cost_ref, u_norm_ref,
1855                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1856
1857        # Check blocking factors with scaling
1858        opts['init_traj'] = res
1859        opts['nominal_traj'] = res
1860        res = op.optimize(self.algorithm, opts)
1861        assert_results(res, cost_ref, u_norm_ref,
1862                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1863
1864        # Check blocking factors for non-inputs
1865        du_bounds = {'cstr.Tc': 15}
1866        factors = {'cstr.Tc': blocking_factors}
1867        bf = BlockingFactors(factors, du_bounds=du_bounds)
1868        opts['blocking_factors'] = bf
1869        N.testing.assert_raises(ValueError, op.optimize, self.algorithm, opts)
1870
1871        # Check blocking factors with bounds and scaling
1872        du_bounds = {'u': 15}
1873        factors = {'u': blocking_factors}
1874        bf = BlockingFactors(factors, du_bounds=du_bounds)
1875        opts['blocking_factors'] = bf
1876        res = op.optimize(self.algorithm, opts)
1877        assert_results(res, 2.280259e3, 3.07723e2,
1878                       cost_rtol=1e-2, u_norm_rtol=5e-3)
1879
1880    @testattr(casadi_base_new = True)
1881    def test_eliminate_der_var(self):
1882        """
1883        Test that results are consistent regardless of eliminate_der_var.
1884        """
1885        mayer_op = self.vdp_bounds_mayer_op
1886        lagrange_op = self.vdp_bounds_lagrange_op
1887       
1888        # References values
1889        cost_ref = 3.17619580332244e0
1890        u_norm_ref = 2.8723837585e-1
1891       
1892        # Keep derivative variables
1893        opts = self.optimize_options(lagrange_op, self.algorithm)
1894        opts["eliminate_der_var"] = False
1895        res = lagrange_op.optimize(self.algorithm, opts)
1896        assert_results(res, cost_ref, u_norm_ref)
1897       
1898        # Mayer, eliminate derivative variables
1899        opts["eliminate_der_var"] = True
1900        opts['init_traj'] = ResultDymolaTextual(
1901                "VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt")
1902        res = mayer_op.optimize(self.algorithm, opts)
1903        assert_results(res, cost_ref, u_norm_ref)
1904       
1905        # Kagrange, eliminate derivative variables
1906        res = lagrange_op.optimize(self.algorithm, opts)
1907        assert_results(res, cost_ref, u_norm_ref)
1908
1909    @testattr(casadi_base_new = True)
1910    def test_eliminate_cont_var(self):
1911        """
1912        Test that results are consistent regardless of eliminate_cont_var.
1913       
1914        This is tested for both Gauss and Radau collocation.
1915        """
1916        op = self._vdp_bounds_mayerop
1917       
1918        # References values
1919        cost_ref = 3.17619580332244e0
1920        u_norm_ref_radau = 2.8723837585e-1
1921        u_norm_ref_gauss = 2.852405405154352e-1
1922       
1923        # Keep continuity variables, Radau
1924        opts = self.optimize_options(op, self.algorithm)
1925        opts['discr'] = "LGR"
1926        opts["eliminate_cont_var"] = False
1927        res = op.optimize(self.algorithm, opts)
1928        assert_results(res, cost_ref, u_norm_ref_radau)
1929       
1930        # Eliminate continuity variables, Radau
1931        opts["eliminate_cont_var"] = True
1932        opts['init_traj'] = ResultDymolaTextual(
1933                "VDP_pack_VDP_Opt_Bounds_Mayer_result.txt")
1934        res = op.optimize(self.algorithm, opts)
1935        assert_results(res, cost_ref, u_norm_ref_radau)
1936       
1937        # Keep continuity variables, Gauss
1938        opts['discr'] = "LG"
1939        opts["eliminate_cont_var"] = False
1940        res = op.optimize(self.algorithm, opts)
1941        assert_results(res, cost_ref, u_norm_ref_gauss)
1942       
1943        # Eliminate continuity variables, Gauss
1944        opts["eliminate_cont_var"] = True
1945        opts['init_traj'] = ResultDymolaTextual(
1946                "VDP_pack_VDP_Opt_Bounds_Mayer_result.txt")
1947        res = op.optimize(self.algorithm, opts)
1948        assert_results(res, cost_ref, u_norm_ref_gauss)
1949
1950    @testattr(casadi_base = True)
1951    def test_quadrature_constraint(self):
1952        """
1953        Test that optimization results of the CSTR is consistent regardless of
1954        quadrature_constraint and eliminate_cont_var for Gauss collocation.
1955        """
1956        op = self.cstr_mayer_op
1957       
1958        # References values
1959        cost_ref = 1.8576873858261e3
1960        u_norm_ref = 3.050971000653911e2
1961       
1962        # Quadrature constraint, with continuity variables
1963        opts = self.optimize_options(op, self.algorithm)
1964        opts['discr'] = "LG"
1965        opts['quadrature_constraint'] = True
1966        opts['eliminate_cont_var'] = False
1967        res = op.optimize(self.algorithm, opts)
1968        assert_results(res, cost_ref, u_norm_ref)
1969       
1970        #~ # Quadrature constraint, without continuity variables #3509
1971        #~ opts['quadrature_constraint'] = True
1972        #~ opts['eliminate_cont_var'] = True
1973        #~ res = op.optimize(self.algorithm, opts)
1974        #~ assert_results(res, cost_ref, u_norm_ref)
1975       
1976        # Evaluation constraint, with continuity variables
1977        opts['quadrature_constraint'] = False
1978        opts['eliminate_cont_var'] = False
1979        res = op.optimize(self.algorithm, opts)
1980        assert_results(res, cost_ref, u_norm_ref)
1981       
1982        #~ # Evaluation constraint, without continuity variables #3509
1983        #~ opts['quadrature_constraint'] = False
1984        #~ opts['eliminate_cont_var'] = True
1985        #~ res = op.optimize(self.algorithm, opts)
1986        #~ assert_results(res, cost_ref, u_norm_ref)
1987
1988    @testattr(casadi_base = True)
1989    def test_n_cp(self):
1990        """
1991        Test varying n_e and n_cp.
1992        """
1993        op = self.vdp_bounds_mayer_op
1994        opts = self.optimize_options(op, self.algorithm)
1995       
1996        # n_cp = 1
1997        opts['n_e'] = 100
1998        opts['n_cp'] = 1
1999        res = op.optimize(self.algorithm, opts)
2000        assert_results(res, 2.58046279958e0, 2.567510746260e-1)
2001       
2002        # n_cp = 3
2003        opts['n_e'] = 50
2004        opts['n_cp'] = 3
2005        res = op.optimize(self.algorithm, opts)
2006        assert_results(res, 3.1761957722665e0, 2.87238440058e-1)
2007       
2008        # n_cp = 8
2009        opts['n_e'] = 20
2010        opts['n_cp'] = 8
2011        opts['init_traj'] = ResultDymolaTextual(
2012                "VDP_pack_VDP_Opt_Bounds_Mayer_result.txt")
2013        res = op.optimize(self.algorithm, opts)
2014        assert_results(res, 3.17620203643878e0, 2.803233013e-1)
2015
2016    @testattr(casadi_base_new = True)
2017    def test_graphs_and_hessian(self):
2018        """
2019        Test that results are consistent regardless of graph and Hessian.
2020       
2021        The test also checks the elimination of derivative and continuity
2022        variables.
2023        """
2024        op = self.vdp_bounds_lagrange_op
2025       
2026        # References values
2027        cost_ref = 3.17619580332244e0
2028        u_norm_ref = 2.8723837585e-1
2029       
2030        # Solve problem to get initialization trajectory
2031        opts = self.optimize_options(op, self.algorithm)
2032        res = op.optimize(self.algorithm, opts)
2033        assert_results(res, cost_ref, u_norm_ref)
2034        opts['init_traj'] = ResultDymolaTextual(
2035                "VDP_pack_VDP_Opt_Bounds_Lagrange_result.txt")
2036       
2037        # SX with exact Hessian and eliminated variables
2038        opts['graph'] = "SX"
2039        opts['IPOPT_options']['hessian_approximation'] = "exact"
2040        opts['eliminate_der_var'] = True
2041        opts['eliminate_cont_var'] = True
2042        res = op.optimize(self.algorithm, opts)
2043        sol_with = res.times['sol']
2044        assert_results(res, cost_ref, u_norm_ref)
2045       
2046        # SX without exact Hessian and eliminated variables
2047        opts['IPOPT_options']['hessian_approximation'] = "limited-memory"
2048        opts['eliminate_der_var'] = False
2049        opts['eliminate_cont_var'] = False
2050        res = op.optimize(self.algorithm, opts)
2051        sol_without = res.times['sol']
2052        nose.tools.assert_true(sol_with < 0.9 * sol_without)
2053        assert_results(res, cost_ref, u_norm_ref)
2054       
2055        # Expanded MX with exact Hessian and eliminated variables
2056        opts['graph'] = "MX"
2057        opts['IPOPT_options']['expand'] = True
2058        opts['IPOPT_options']['hessian_approximation'] = "exact"
2059        opts['eliminate_der_var'] = True
2060        opts['eliminate_cont_var'] = True
2061        res = op.optimize(self.algorithm, opts)
2062        sol_with = res.times['sol']
2063        assert_results(res, cost_ref, u_norm_ref)
2064       
2065        # Expanded MX without exact Hessian and eliminated variables
2066        opts['IPOPT_options']['hessian_approximation'] = "limited-memory"
2067        opts['eliminate_der_var'] = False
2068        opts['eliminate_cont_var'] = False
2069        res = op.optimize(self.algorithm, opts)
2070        sol_without = res.times['sol']
2071        nose.tools.assert_true(sol_with < 0.9 * sol_without)
2072        assert_results(res, cost_ref, u_norm_ref)
2073       
2074        # MX with exact Hessian and eliminated variables
2075        opts['IPOPT_options']['expand'] = False
2076        opts['IPOPT_options']['hessian_approximation'] = "exact"
2077        opts['eliminate_der_var'] = True
2078        opts['eliminate_cont_var'] = True
2079        res = op.optimize(self.algorithm, opts)
2080        assert_results(res, cost_ref, u_norm_ref)
2081       
2082        # MX without exact Hessian and eliminated variables
2083        opts['IPOPT_options']['hessian_approximation'] = "exact"
2084        opts['eliminate_der_var'] = False
2085        opts['eliminate_cont_var'] = False
2086        res = op.optimize(self.algorithm, opts)
2087        assert_results(res, cost_ref, u_norm_ref)
2088
2089    @testattr(casadi_base = True)
2090    def test_solver_statistics(self):
2091        """
2092        Test getting IPOPT statistics
2093        """
2094        op = self.vdp_bounds_mayer_op
2095       
2096        cost_ref = 3.17619580332244e0
2097       
2098        res = op.optimize()
2099        (return_status, nbr_iter, objective, total_exec_time) = \
2100                res.get_solver_statistics()
2101        N.testing.assert_string_equal(return_status, "Solve_Succeeded")
2102        N.testing.assert_array_less([nbr_iter, -nbr_iter], [100, -5])
2103        N.testing.assert_allclose(objective, cost_ref, 1e-3, 1e-4)
2104        N.testing.assert_array_less([total_exec_time, -total_exec_time],
2105                                    [10., 0.])
2106
2107    @testattr(casadi_base = True)
2108    def test_input_interpolator(self):
2109        """
2110        Test the input interpolator for simulation purposes
2111        """
2112        model = self.cstr_model
2113        model.reset()
2114        op = self.cstr_extends_op
2115        model.set(['c_init', 'T_init'], op.get(['c_init', 'T_init']))
2116       
2117        # Optimize
2118        opts = self.optimize_options(op, self.algorithm)
2119        opts['n_e'] = 100       
2120        opt_res = op.optimize(self.algorithm, opts)
2121       
2122        # Simulate
2123        opt_input = opt_res.get_opt_input()
2124        opts = model.simulate_options()
2125        opts["CVode_options"]["rtol"] = 1e-6
2126        opts["CVode_options"]["atol"] = 1e-8 * model.nominal_continuous_states
2127        res = model.simulate(start_time=0., final_time=150., input=opt_input,
2128                             options=opts)
2129        N.testing.assert_allclose([res.final("T"), res.final("c")],
2130                                  [284.60202203, 346.31140851], rtol=5e-4)
2131
2132    @testattr(casadi_base = True)
2133    def test_matrix_evaluations(self):
2134        """
2135        Test evaluating NLP matrices.
2136        """
2137        op = self.cstr_lagrange_op
2138
2139        # References values
2140        cost_ref = 1.852527678e3
2141        u_norm_ref = 3.045679e2
2142
2143        # Solve
2144        opts = self.optimize_options(op, self.algorithm)
2145        opts['n_e'] = 20
2146        res = op.optimize(self.algorithm, opts)
2147        assert_results(res, cost_ref, u_norm_ref, u_norm_rtol=5e-3)
2148
2149        # Compute Jacobian condition numbers
2150        J_init = res.solver.get_J("init")
2151        J_init_cond = N.linalg.cond(J_init)
2152        N.testing.assert_allclose(J_init_cond, 3.3368e5, rtol=1e-2) #3.3368e5 , #2.93e4
2153        J_opt = res.solver.get_J("opt")
2154        J_opt_cond = N.linalg.cond(J_opt)
2155        N.testing.assert_allclose(J_opt_cond, 3.6453652e7, rtol=1e-2) #3.6453652e7 , #3.37e6
2156
2157        # Compute Hessian norms
2158        H_init = res.solver.get_H("init")
2159        H_init_norm = N.linalg.norm(H_init)
2160        N.testing.assert_allclose(H_init_norm, 4.36e3, rtol=1e-2)
2161        H_opt = res.solver.get_H("opt")
2162        H_opt_norm = N.linalg.norm(H_opt)
2163        N.testing.assert_allclose(H_opt_norm, 1.47e5, rtol=1e-2)
2164
2165        # Compute KKT condition numbers
2166        KKT_init = res.solver.get_KKT("init")
2167        KKT_init_cond = N.linalg.cond(KKT_init)
2168        N.testing.assert_allclose(KKT_init_cond, 2.637e9, rtol=1e-2) #2.637e9 #2.72e8
2169        KKT_opt = res.solver.get_KKT("opt")
2170        KKT_opt_cond = N.linalg.cond(KKT_opt)
2171        N.testing.assert_allclose(KKT_opt_cond, 9.705e9, rtol=1e-2) #9.705e11 #1.18e10
2172
2173        # Obtain symbolic matrices and matrix functions
2174        res.solver.get_J("sym")
2175        res.solver.get_J("fcn")
2176        res.solver.get_H("sym")
2177        res.solver.get_H("fcn")
2178        res.solver.get_KKT("sym")
2179        res.solver.get_KKT("fcn")
2180
2181    @testattr(casadi_base = True)
2182    def test_expand_to_sx(self):
2183        """
2184        Test the different levels of SX expansion.
2185        """
2186        op = self.vdp_bounds_lagrange_op
2187
2188        # Reference values
2189        cost_ref = 3.17619580332244e0
2190        u_norm_ref = 2.8723837585e-1
2191        opts = self.optimize_options(op)
2192
2193        # Test NLP expansion
2194        opts['expand_to_sx'] = "NLP"
2195        res = op.optimize(options=opts)
2196        assert_results(res, cost_ref, u_norm_ref)
2197       
2198        # Test DAE expansion
2199        opts['expand_to_sx'] = "DAE"
2200        res = op.optimize(options=opts)
2201        assert_results(res, cost_ref, u_norm_ref)
2202
2203        # Test without expansion
2204        opts['expand_to_sx'] = "no"
2205        res = op.optimize(options=opts)
2206        assert_results(res, cost_ref, u_norm_ref)
2207
2208    @testattr(casadi_base = True)
2209    def test_vdp_function(self):
2210        """
2211        Test a VDP model containing a function.
2212        """
2213        op = self.vdp_function_op
2214
2215        # Reference values
2216        cost_ref = 3.17619580332244e0
2217        u_norm_ref = 2.8723837585e-1
2218
2219        opts = self.optimize_options(op)
2220
2221        # Test with full SX expansion
2222        opts['expand_to_sx'] = "NLP"
2223        res = op.optimize(options=opts)
2224        assert_results(res, cost_ref, u_norm_ref)
2225
2226        # Test without SX expansion
2227        opts['expand_to_sx'] = "no"
2228        res = op.optimize(options=opts)
2229        assert_results(res, cost_ref, u_norm_ref)
2230
2231    @testattr(casadi_base = True)
2232    def test_vdp_function_global_constant(self):
2233        """
2234        Test a VDP model containing a function with global constant.
2235        """
2236        op = self.vdp_function_op_global_constant
2237
2238        # Reference values
2239        cost_ref = 3.17619580332244e0
2240        u_norm_ref = 2.8723837585e-1
2241
2242        opts = self.optimize_options(op)
2243
2244        # Test with full SX expansion
2245        opts['expand_to_sx'] = "NLP"
2246        res = op.optimize(options=opts)
2247        assert_results(res, cost_ref, u_norm_ref)
2248
2249        # Test without SX expansion
2250        opts['expand_to_sx'] = "no"
2251        res = op.optimize(options=opts)
2252        assert_results(res, cost_ref, u_norm_ref)
2253
2254    @testattr(worhp = True)
2255    def test_worhp(self):
2256        """
2257        Test using WORHP instead of IPOPT.
2258        """
2259        op = self.vdp_bounds_lagrange_op
2260        opts = self.optimize_options(op)
2261
2262        # Reference values
2263        cost_ref = 3.17619580332244e0
2264        u_norm_ref = 2.8723837585e-1
2265
2266        # Test with IPOPT
2267        opts['solver'] = 'IPOPT'
2268        res = op.optimize(options=opts)
2269        assert_results(res, cost_ref, u_norm_ref)
2270
2271        # Test with WORHP
2272        opts['solver'] = 'WORHP'
2273        opts['WORHP_options']['NLPprint'] = -1
2274        res = op.optimize(options=opts)
2275        assert_results(res, cost_ref, u_norm_ref)
2276   
2277    @testattr(casadi_base = True)
2278    def test_instream(self):
2279        """
2280        Test inStream operator support.
2281        """
2282        op = self.stream_op
2283        opts = self.optimize_options(op)
2284       
2285        # Set arbitrary non-zero guesses
2286        guesses = N.linspace(0.5, 1.8, len(op.getAllVariables()))
2287        for (var, ig) in zip(op.getAllVariables(), guesses):
2288            var.setAttribute('initialGuess', ig)
2289
2290        # Reference values
2291        cost_ref = 63.99032036914
2292        u_norm_ref = 5.9669830727
2293       
2294        # Solve and verify result
2295        opts['solver'] = 'IPOPT'
2296        opts['n_e'] = 3
2297        opts['named_vars'] = True
2298        res = op.optimize(options=opts)
2299        assert_results(res, cost_ref, u_norm_ref)
2300
2301class TestLocalDAECollocator_expand_to_sx_DAE(TestLocalDAECollocator):
2302    """
2303    Tests pyjmi.optimization.casadi_collocation.LocalDAECollocator with option 'expand_to_sx': 'DAE'
2304    """
2305    def optimize_options(self, *args, **kwargs):
2306        opts = TestLocalDAECollocator.optimize_options(self, *args, **kwargs)
2307        opts['expand_to_sx'] = 'DAE'
2308        return opts
2309
2310    @nose.tools.nottest
2311    def test_cstr_checkpoint(self):
2312        assert False # shouldn't be run, just used to disable this test with 'expand_to_sx': 'DAE'
2313
2314class TestLocalDAECollocator_BLT_transfer(TestLocalDAECollocator):
2315    """
2316    Tests pyjmi.optimization.casadi_collocation.LocalDAECollocator with option 'equation_sorting':True
2317    """
2318   
2319    @staticmethod
2320    def compile_options():
2321        compiler_options = {'equation_sorting':True}
2322        return compiler_options
2323   
2324    @nose.tools.nottest
2325    def test_cstr_checkpoint(self):
2326        assert False # shouldn't be run, just used to disable this test with 'expand_to_sx': 'DAE'   
2327   
2328    @nose.tools.nottest
2329    def test_result_modes(self):
2330        assert False
2331   
2332    @nose.tools.nottest   
2333    def test_cstr_minimum_time(self):
2334        assert False
2335       
2336    @nose.tools.nottest   
2337    def test_result_mode(self):
2338        assert False 
2339       
2340    @nose.tools.nottest   
2341    def test_vdp_minimum_time(self):
2342        assert False
2343       
2344    @nose.tools.nottest   
2345    def test_quadrature_constraint(self):
2346        assert False 
2347       
2348    @nose.tools.nottest   
2349    def test_parameter_estimation_traj(self):
2350        assert False   
2351       
2352    @nose.tools.nottest   
2353    def test_parameter_estimation(self):
2354        assert False 
2355   
2356    @nose.tools.nottest   
2357    def test_parameter_setting(self):
2358        assert False 
2359       
2360    @nose.tools.nottest   
2361    def test_nominal_traj_mode(self):
2362        assert False 
2363       
2364    @nose.tools.nottest   
2365    def test_matrix_evaluations(self):
2366        assert False 
2367       
2368    @nose.tools.nottest   
2369    def test_nominal_traj_cstr(self):
2370        assert False 
2371       
2372    @nose.tools.nottest   
2373    def test_nominal_traj_vdp(self):
2374        assert False 
2375
2376    def test_named_vars_c_e_10_id(self):
2377        return 11;
Note: See TracBrowser for help on using the repository browser.