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

Last change on this file since 13800 was 13800, checked in by randersson, 8 weeks ago

#5819 Merged trunk into branch

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