SÃO PAULO STATE UNIVERSITY “JÚLIO DE MESQUITA FILHO” CAMPUS BAURU SCHOOL OF ENGINEERING RANULFO ACIR DE OLIVEIRA RESENDE Nonconvexity in Optimal Control Problems: An Approach by Automated Numerical Simulations Bauru - SP 2018 RANULFO ACIR DE OLIVEIRA RESENDE Nonconvexity in Optimal Control Problems: An Approach by Automated Numerical Simulations Thesis presented to São Paulo State University (UNESP), Campus Bauru, School of Engineering, as part of the requirements for obtaining the title of Ph.D. in ELECTRICAL ENGINEERING. Prof. Tit. José Manoel Balthazar Advisor Bauru - SP 2018 Resende, Ranulfo Acir de Oliveira. Nonconvexity in optimal control problems : an approach by automated numerical simulations / Ranulfo Acir de Oliveira Resende, 2018 142 f. : il. Advisor: José Manoel Balthazar Thesis (Doctoral)–São Paulo State University (UNESP). School of Engineering, Bauru, 2018 1. Optimization. 2. Optimal control. 3. Mathematical modeling. 4. Numerical simulation. I. São Paulo State University (UNESP). School of Engineering. II. Title. Acknowledgements My gratitude: • To God, for giving me the strength and health to carry out this work; • To my family, for the support and encouragement; • To Prof. Tit. José Manoel Balthazar, for the academic guidance; and • To the School of Engineering, Campus Bauru, Unesp, its faculty, staff and administration, for the support to this research. “O mundo é um moinho.” Agenor de Oliveira - Cartola (1908-1980) “No mundo passais por aflições; mas tende bom ânimo, eu venci o mundo.” Bíblia (João 16:33) À minha família, de quem aprendemos e pra quem ensinamos o amor, o perdão, a esperança e a fé, com quem aprendemos e ensinamos o amor, o perdão, a esperança e a fé. Resumo Esta Tese de Doutorado é composta pelo estudo da aplicação de simulações numéricas automatizadas para solucionar problemas de controle ótimo (PCO) não convexos. Um algoritmo para tratar PCOs não convexos é apresentado como a principal contribuição da Tese. Outra contribuição, a estrutura ASAF - Automated Simulation and Analysis Framework, codificada em linguagem Python e baseada em modernas ferramentas de modelagem, simulação, JModelica.org, e otimização, Interior-Point Optimizer (IPOPT), é utilizada como auxílio para aplicação do referido algorítmo. O algoritmo proposto foi testado de forma metodológica com escala métrica e quatro PCOs não convexos com soluções analíticas exatas conhecidas. A viabilidade de aplicação do algoritmo foi verificada por clássicos métodos estatísticos aplicados aos dados das simulações, alcançando-se os resultados esperados e fundamentando-se a possibilidade de se tratar PCOs não convexos de soluções desconhecidas. Assim, esta Tese apresenta uma contribuição relevante para o tratamento de PCOs não convexos, ao propor um algoritmo para, com uso eficaz dos crescentes recursos computacionais de metaprogramação e automatização, obter solução subótima recomendada aceitável. A mediana da efetividade do algoritmo é estimada, a partir de testes estatísticos, em 93%, com significância de 0,05 (95% de confiança). Palavras-chave: otimização. controle ótimo. modelagem matemática. simulação numérica. Abstract This doctoral Thesis is comprised by the study of the application of automated numerical simulations to solve nonconvex optimal control problems (OCP). An algorithm to treat nonconvex OCPs is presented as the main contribution of the Thesis. Another contribution of the Thesis, the ASAF-Automated Simulation and Analysis Framework, coded in Python language and based on modern tools of modeling, simulation, JModelica.org, and optimization, Interior-Point Optimizer (IPOPT), is used as an aid in the application of this algorithm. The proposed algorithm was tested by a methodological procedure with a metric scale and four nonconvex OCPs with known exact analytical solutions. The feasibility of using the algorithm was verified by classic statistical methods applied to the simulations data, achieving the expected outcome and grounding the possibility to treat nonconvex OCPs of unknown solutions. Thus, this Thesis presents a relevant contribution to the treatment of nonconvex OCPs, by proposing an algorithm to, upon the effective use of the increasing computer resources of metaprogramming and automation, obtain an acceptable sub-optimal recommended solution. The median of the algorithm effectiveness is estimated, from statistical tests, in 93%, with a significance of 0.05 (95% confidence). Keywords: Optimization. Optimal control. Mathematical modeling. Numerical simulation. List of Figures Figure 1 – Graphic Description of the Contribution of the Thesis for the Treatment of Nonconvex OCPs. Source: The author. 27 Figure 2 – Sequence for the Application of the Algorithm. Source: The author. 44 Figure 3 – Description of the Process of Algorithm Application. Source: The author. 48 Figure 4 – Automated Structure to Solve Nonlinear and Nonconvex OCPs. Source: The author. 50 Figure 5 – Development of the Graphical User Interface. Source: The author 56 Figure 6 – First Analytical Solution of OCP#1. Source: The author. 60 Figure 7 – Second Analytical Solution of OCP#1. Source: The author. 61 Figure 8 – First Analytical Solution for OCP#2. Source: The author. 63 Figure 9 – Second Analytical Solution for OCP#2. Source: The author. 63 Figure 10 – Candidate Solution for OCP#1 - Context A: Control. Source: The author. 68 Figure 11 – Candidate Solution for OCP#1 - Context A: State. Source: The author. 69 Figure 12 – Rejected Solution for OCP#1 - Context A: Control. Source: The author. 70 Figure 13 – Rejected Solution for OCP#1 - Context A: State. Source: The author. 70 Figure 14 – Solution Improperly Classified as Candidate for OCP#1 - Context A: Control. Source: The author. 71 Figure 15 – Solution Improperly Classified as Candidate for OCP#1- Context A: State. Source: The author. 72 Figure 16 – Candidate Solution for OCP#1 - Context B: Control. Source: The author. 73 Figure 17 – Candidate Solution for OCP#1 - Context B: State. Source: The author. 74 Figure 18 – Rejected Solution for OCP#1 - Context B: Control. Source: The author. 75 Figure 19 – Rejected Solution for OCP#1 - Context B: State. Source: The author. 75 Figure 20 – Candidate Solution for OCP#1 - Context C: Control. Source: The author. 76 Figure 21 – Candidate Solution for OCP#1 - Context C: State. Source: The author. 76 Figure 22 – Candidate Solution for OCP#1 - Context S: Control. Source: The author. 77 Figure 23 – Candidate Solution for OCP#1 - Context S: State. Source: The author. 78 Figure 24 – Rejected Solution for OCP#1 - Context S: Control. Source: The author. 79 Figure 25 – Rejected Solution for OCP#1 - Context S: Control. Source: The author. 79 Figure 26 – Candidate Solution for OCP#1 - Context T: Control. Source: The author. 80 Figure 27 – Candidate Solution for OCP#1 - Context T: State. Source: The author. 80 Figure 28 – Rejected Solution for OCP#1 - Context T: Control. Source: The author. 81 Figure 29 – Rejected Solution for OCP#1 - Context T: State. Source: The author. 81 Figure 30 – Candidate Solution for OCP#1 - Context U: Control. Source: The author. 82 Figure 31 – Candidate Solution for OCP#1 - Context U: State. Source: The author. 82 Figure 32 – Rejected Solution for OCP#1 - Context U: Control. Source: The author. 83 Figure 33 – Rejected Solution for OCP#1 - Context U: State. Source: The author. 83 List of Tables Table 1 – Results of Searches in Databases 33 Table 2 – Searches in Databases: Most Cited Articles on the Topic 34 Table 3 – Parameters for the Solution Classification as Candidate 43 Table 4 – Tools Used in this Thesis. 52 Table 5 – Summary of the Canonical Solutions of the Used OCPs 65 Table 6 – Optimization Parameters Used in All the Optimizations 65 Table 7 – Parameters Used in the Definitions of the Contexts 66 Table 8 – Parameters Used in the Definitions of the Simulations 67 Table 9 – Scope of Simulations 67 Table 10 – Optimizations Results for OCP#1 67 Table 11 – Estimation of the Effectiveness of the Algorithm: Summary of the Results 84 List of Abbreviations and Acronyms HJB Hamilton-Jacobi-Bellman (Equation of) IPOPT Interior-Point Optimizer (Optimization Tool) MPC Model Predictive Control NLP Nonlinear Program NMPC Nonlinear Model Predictive Control NRMSE Normalized Root-Mean-Squared Error OCP Optimal Control Problem OP Optimization Problem PMP Pontryagin Maximum Principle List of Symbols α Significance level1 α Maximum non-compliance index 2,3 ρ Maximum relative rate for growth of J 2 ε Minimum acceptable slack to the bounds of the variables2 uuu0 Initial condition for control variable uuu∗ Optimal control solution uuu Control variable - control vector xxx0 Initial condition for state variable xxx∗ Optimal state solution xxx State variable - state vector t0 Initial condition for time variable t f Final condition for time variable J[xxx,uuu, t] Performance Index, Cost Function or Objective ne Number of elements for the discretization of the OCP4 nc Number of collocation points in each element4 τi,k Collocation point ’k’ in element ’i’4 ˜̀k Lagrange basis polynomials with interpolation point τi,0 4 `k Lagrange basis polynomials without interpolation point τi,0 4 gi Inequality constraints4 ge Equality constraints4 Gi Point inequality constraints4 Ge Point equality constraints4 1Statistical Parameter. 2Parameter for the Analysis of the Numerical Solution. 3Indicates the nonconvexity of solution uuu∗ or xxx∗, non-monotonicity or non-slow growth of J. 4Applicable in the technique of dynamic optimization by collocation points. Contents 1 INTRODUCTION 25 2 CONSIDERATIONS ABOUT THE METHODOLOGY 31 2.1 BIBLIOGRAPHICAL REVIEW 31 2.2 THE SCIENTIFIC PROBLEM 35 2.3 THE HYPOTHESIS 35 2.4 THE REFUTATION METHOD 36 2.5 STATISTICAL ANALYSIS AND TOOLS TO VERIFY THE HYPOTHESIS 38 3 CONTRIBUTIONS 41 3.1 ALGORITHM TO TREAT NONCONVEX OPTIMAL CONTROL PROBLEMS 41 3.2 THE INTEGRATED AND AUTOMATED ENVIRONMENT FOR MODELING, SIMULATION AND OPTIMIZATION 48 3.2.1 JModelica.org 49 3.2.2 IPOPT 50 3.2.3 HSL Routines 51 3.2.4 Python and other tools 51 3.3 ASAF - AUTOMATED SIMULATION AND ANALYSIS FRAMEWORK 52 3.3.1 Main Features of the ASAF Structure 54 4 THE HYPOTHESIS VERIFICATION 59 4.1 FIRST OCP USED IN THE HYPOTHESIS VERIFICATION 59 4.2 SECOND OCP USED IN THE HYPOTHESIS VERIFICATION 62 4.3 THIRD OCP USED IN THE HYPOTHESIS VERIFICATION 64 4.4 FOURTH OCP USED IN THE HYPOTHESIS VERIFICATION 64 4.5 RESULTS 65 4.5.1 Optimizations Results for OCP#1 67 4.5.1.1 Analysis of the Optimizations Results for OCP#1 68 4.6 THE HYPOTHESIS VERIFICATION BY STATISTICAL TESTS 84 4.6.1 The MATLAB R© Code for the Hypothesis Verification 84 5 DISCUSSION OF THE RESULTS 87 6 CONCLUSION 91 REFERENCES 93 APPENDIX A SUMMARY OF THE DYNAMIC OPTIMIZATION BY DIRECT COLLOCATION 99 A.1 COLLOCATION POLYNOMIALS 102 A.2 THE TRANSFORMATION OF THE OCP 103 APPENDIX B STATISTICAL ANALYSIS OF THE SIMULATIONS 107 APPENDIX C THE OUTPUT OF THE KRUSKAL-WALLIS TEST 135 INDEX 141 25 1 INTRODUCTION The theory of optimal control aims to minimize a performance index which, most of the time, is related to a function of economic cost. This theory is the core of this research and there are many excellent bibliographical sources, such as (KIRK, 1970; BRYSON, 1975; ATHANS, 2007), with no intention of discrediting others, that present various methods and techniques to solve specific categories of Optimal Control Problems (OCPs). Classical techniques and methods to solve OCPs rely on the Theorems of Existence (CESARI, 1983, p. 367-472) to ensure optimal solution. However, there is a variety of OCPs that do not comply with the classical conditions of linearity and convexity. Currently, the theory of optimal control shows relevant evolution in the theorems of existence, for which the non-compliance with the classical conditions (CESARI, 1974), (ANCONA, 2008), (CELLINA; COLOMBO, 1990), (MARICONDA, 1992), (EKELAND; TEMAM, 1976), (ZASLAVSKI, 2007), (CLARKE, 1983), (CLARKE, 2011), such as convexity and lower semi-continuity (CLARKE, 1983, p.103), under certain conditions, does not compromise the existence of optimal solution, even if it is very hard to find it. With the emergence of recent concepts and approaches defined by new conditions of existence (ZASLAVSKI, 2013, p. 1) based on the performance index growth, the theoretical evolution was not followed by the practical development to obtain optimal solutions. This Thesis seeks to fill this gap, looking to use the potential of metaprogramming, 26 Chapter 1. INTRODUCTION the versatility and flexibility of the models in Modelica language and computer automation in numerical simulations, to advance the study of nonconvex OCPs and exploit these new concepts. 27 Pr ob le m s T re at ed b y N um er ic al S im ul at io ns C on tr ol P ro bl em s O C P s C om pl ia nt O C P s N on -C om pl ia nt O C P s Fi gu re 1 – G ra ph ic D es cr ip tio n of th e C on tr ib ut io n of th e T he si s fo rt he Tr ea tm en to fN on co nv ex O C Ps .S ou rc e: T he au th or . 28 Chapter 1. INTRODUCTION Figure 1 graphically describes the contribution of this Thesis to the expansion of the nonconvex set of OCPs with acceptable solutions obtained by numerical simulations. Therefore, the overall goal of this research was to find a method to treat non-compliant OCPs and obtain acceptable solutions through automated numerical simulations. This research presents as main result, in Section 3.1, the definition of an algorithm to achieve acceptable solutions for non-compliant OCPs (nonconvex, nonlinear) using the assistance of automated numerical simulations. In Section 3.3, the ASAF - Automated Simulation and Analysis Framework, a set of Python scripts used to verify the suggested algorithm is presented as secondary result. The median of the algorithm effectiveness was estimated, from statistical tests, in 93%, with a significance of 0.05 (95% confidence). The process to verify the effectiveness of the algorithm carried out 972 simulations of four OCPs. These OCPs were specifically chosen from the literature because they present nonconvexity and known solutions for state vector. Therefore, a relevant contribution to the treatment of nonconvex OCPs is provided by this Thesis. The Thesis was divided into the following chapters. Chapter 2, CONSIDERATIONS ABOUT THE METHODOLOGY, presents the relevant topics that were followed to achieve the Thesis results: Initially, the details of the bibliographical and theoretical review are described, and, after that, the Scientific Problem, the Hypothesis and The Refutation Method are defined. Then, in Chapter 3, CONTRIBUTIONS, the main and the secondary results of the Thesis are presented, along with a short description of the software tools used. Chapter 4, THE HYPOTHESIS VERIFICATION, contains the details about the procedures performed to verify the hypothesis of this Thesis. Chapter 5, DISCUSSION OF THE RESULTS, discusses the implications of the achieved 29 results in relation to the current knowledge about the treatment of non-compliant OCPs, as well as to the possibilities of evolution of computer resources and of automated numerical simulations. Chapter 6, CONCLUSION, presents a brief description of the importance of the achieved results. Next, in REFERENCES, all bibliographical publications used to support this doctoral Thesis are provided. The Appendix A presents a SUMMARY OF THE DYNAMIC OPTIMIZATION BY DIRECT COLLOCATION, used in the automated simulations, based on the original source (MAGNUSSON; ÅKESSON, 2012). A report of the numerical results and the solutions classifications is presented in Appendix B. Appendix C presents the statistical analysis of the Kruskal-Wallis Test (DANIEL, 1990), which estimated the median of the algorithm effectiveness. The next chapter presents the details about the methodology used in this Ph.D. Thesis. 31 2 CONSIDERATIONS ABOUT THE METHODOLOGY This chapter presents the methodological details, (MARCONI; LAKATOS, 2010; KOCHE, 2009; ECO, 2007), that have guided the development of the Thesis and explains how the results were achieved. Initially, the procedures performed within the bibliographical review are described. Next, the Scientific Problem, the Hypothesis and the Refutation Method used in the Thesis are presented. 2.1 BIBLIOGRAPHICAL REVIEW A bibliographical review was carried out, which started with classical literature on Control Engineering, (OGATA, 2010; OGATA, 1995; FRANKLIN; POWELL, 1980; FRANKLIN, 1991), in which some fundamental concepts related essentially to the Modern Control Theory are evidenced, such as state variables and state equations. Next, the optimal control theory was studied, (LEWIS, 1986), (BRYSON, 1975), (KIRK, 1970), (ATHANS, 2007), (BERTSEKAS, 2005a; BERTSEKAS, 2005b), where the performance index, the Quadratic Continuous Linear Regulator, the Linear Quadratic Continuous Tracker, the Pontryagin Maximum Principle, the Bellman Optimality Principle, the Hamilton-Jacobi-Bellman and Riccati equations were defined. Similarly, a bibliographical review and theoretical foundation related to Optimization were performed, as a science that seeks to solve a static optimization problem (OP), may it be 32 Chapter 2. CONSIDERATIONS ABOUT THE METHODOLOGY linear or nonlinear. Classic literature in this area can be found on (LUENBERGER; YE, 2008), (BAZARAA, 2009), where fundamental concepts are presented, such as the Constrained Optimization Problem, Necessary Conditions of First and Secord Order, Karush, Kuhn and Tucker (KKT) Criteria, Optimization Methods, Barrier Method, Interior-Points Method, Lagrange Multipliers, Stopping Criteria, Feasible Region, Feasible Solution, Local and Global optimal solutions. This Thesis used the DYNAMIC OPTIMIZATION BY DIRECT COLLOCATION, which was proposed by (MAGNUSSON; ÅKESSON, 2012), as support for the implementation of automated numerical simulations for solving non-compliant OCPs. The method is summarized in Appendix A. The review proceeded with the use of more recent sources which presented a discussion on nonconvexity, (CESARI, 1974), (ANCONA, 2008), (CLARKE, 1983), (CLARKE, 2011), (CELLINA; COLOMBO, 1990), (MARICONDA, 1992), (EKELAND; TEMAM, 1976) and (ZASLAVSKI, 2007). Next, Existence Theorems (CESARI, 1983, p. 367-472) to ensure the optimal solution were studied. Some of them consider slow and monotonic growth of the performance index and convexity to classify a OCP as solvable. This review has highlighted a relevant fact that, for an eventual non-compliant OCP, for example, nonconvex, the theoretical existence of optimal solution can still be mathematically defined. Zaslaviski shows in (ZASLAVSKI, 2013, p. 2) the existence of a solution for OCP set under the unique condition of characteristics of the objective function with monotonic non-decreasing growth (CESARI, 1983, p. 104). After the cited initial reviews, a survey was performed in some databases available on the Internet, searching for articles on the topics of the Thesis. Searches were carried out in Scopus and Web of Science from the Institute for Scientific 2.1. Bibliographical Review 33 Information (ISI)1. The first search was performed using the keywords "Control Problems" and "Numerical Simulations". In the second search, the expressions "Optimal Control" and "Numerical Simulations" were used. In both cases, the results were filtered for the fields: "Computer Science", "Mathematics" and “Engineering”. The results are presented in Table 1. Table 1 – Results of Searches in Databases Search Scopus ISI Web of Science 1 4034 125 2 4034 536 Source: The author. After inspection of the titles and abstracts of the papers obtained in the searches, with most of them standing out of the scope of this Thesis, a detailed study of the most cited articles was made, with a minimum of 20 citations. The articles were sorted by topics and type of contribution, considering if they were about a theoretical development or an application. A summary of this analysis is presented in the Table 2. 1At https://www.webofknowledge.com/, accessed on may, 6th, 2016. 34 Chapter 2. CONSIDERATIONS ABOUT THE METHODOLOGY Table 2 – Searches in Databases: Most Cited Articles on the Topic # Article Citations Remarks 1 (TODOROV, 2005) 127 development of theory stochastic optimal control coordinate-descendants algorithm 2 (RAFIKOV; BALTHAZAR, 2008) 121 development of theory nonlinear systems chaos synchronization linear feedback controller 3 (LUO, 2005) 106 application tracking control problem 4 (CHESI, 2003) 94 development of theory 5 (JOO, 1999) 91 development of theory digital control of chaotic systems model-based fuzzy controller 6 (SHAO, 2007) 78 application 7 (LIN, 1999) 77 application 8 (BANI-HANI; GHABOUSSI, 1998) 69 development of theory neural networks nonlinear control 9 (BEWLEY, 2000) 66 application linear and nonlinear problems 10 (LI, 2007) 60 development of theory 11 (HOUSKA, 2011) 59 application NMPCa 12 (ILZHOEFER, 2007) 59 application NMPCa 13 (CIMEN; BANKS, 2004) 57 development of theory ASREb 14 (LIU, 2012) 56 application MASc 15 (CHOI, 2004) 53 application 16 (ZHOU, 1995) 47 application 16 (ERKUS, 2002) 45 application 18 (DING, 2004) 40 application 18 (BONNARD, 2009) 39 application 20 (LENCI; REGA, 2003) 35 development of theory a Nonlinear Model Predictive Control b Approximating Sequence Riccati Equations c Multi-agent Systems Source: Survey data. 2.2. The Scientific Problem 35 After this study, it was found that article (RAFIKOV; BALTHAZAR, 2008) is the work which is more closely related to the scope of this Thesis, presenting an approach to the treatment of nonlinearities with linear feedback control. It was also evidenced that there is a gap about the treatment of the problem that is the focus of the Thesis, which means, not compliant with the conditions and classical methods. Another important point about the bibliographical review includes the characteristics of the solutions obtained by numerical simulations. Thus, concepts such as integration of trajectories, integration steps, single-step algorithm, multi-step algorithm, rounding error, truncation error and numerical stability are presented in Chapter 4 of (PARKER; CHUA, 1991). 2.2 THE SCIENTIFIC PROBLEM The scientific problem can be summarized within the following research question: How to treat, through automated numerical simulations, and obtain an acceptable solution for Optimal Control Problems (OCPs) that, despite being interesting, have characteristics of non-compliance (e.g., nonconvexity, nonlinearity) with the classical theory and with the current main solution methods? An algorithm to treat non-compliant OCPs was proposed in Section 3.1, as an answer to the research question 2.3 THE HYPOTHESIS The hypothesis of this doctoral research is defined by: Null hypothesis (H0): The algorithm to treat non-compliant OCPs presented in Section 3.1, when applied in a particular context in which the OCP studied is consistently defined in relation to the theory and to the restrictions of the environment for numeric simulations and optimizations, is able to regularly present a median of effectiveness of correction of solutions classifications equal to or greater than 70%. The alternative hypothesis, that will refute the Thesis, is defined by: 36 Chapter 2. CONSIDERATIONS ABOUT THE METHODOLOGY Alternative hypothesis (H1): The algorithm to treat non-compliant OCPs presented in Section 3.1, when applied in a particular context in which the OCP studied is consistently defined in relation to the theory and to the restrictions of the environment for numeric simulations and optimizations, regularly presents a median of effectiveness of correction of solutions classifications lower than 70%. 2.4 THE REFUTATION METHOD The problem to evaluate and compare solutions from several numerical simulations, in face of the diversity of possible integration algorithms and simulation parameters, is about how to obtain the reference solution of the problem to be compared with the provided solutions. If the OCP in question allows an analytical or exact solution, this solution will be the reference for the comparison. This was the approach found in (PARKER; CHUA, 1991, p. 83), in which various solutions are provided by several integration methods applied to the same problem of a known exact solution. Therefore, in (PARKER; CHUA, 1991, p. 83), the algorithms assessment took place through the comparison of its numerical solution with the exact solution. Obviously, the algorithm proposed in this Thesis does not restrict the OCP in relation to the condition of having a known analytical solution, neither sets restrictions to the objective function. On the contrary, as any other algorithm to be applied with the assistance of numerical simulations, the possible attainment of solutions to generic OCPs with no additional restrictions is precisely its applicability advantage. However, in the context of corroboration of the hypothesis and the obligation to evaluate the solutions provided by the algorithm, one can resort to, as an analogy to that was carried out in (PARKER; CHUA, 1991, p. 83), a comparison technique with the analytical or exact solution. Another consideration is the randomness introduced to the process by variation of parameters performed in steps 2 and 3 of the algorithm, as well as by the diversity of possibilities of parametrically modeling of a given OCP (step 1). 2.4. The Refutation Method 37 Thus, it is considered to be not practical to develop an algorithm with effectiveness of 100% correct classifications with 100% confidence. Therefore, it becomes necessary to elaborate a refutation method that, based on metric indexes (such as Normalized Root-Mean-Squared Error - NRMSE), makes use of statistical methods to test the hypothesis and estimate the median of the effectiveness distribution of the algorithm. The following steps are suggested as the refutation method: 1. Choose an OCP2 with characteristics of non-compliance and that has, in any possible way, an exact known solution; 2. When applying the algorithm presented in Section 3.1, and without loss of generality, associate each region defined on step 4 of the proposed algorithm, chosen by the user based on his or her requirements, to a single context of simulations, as Definition 3.4; 3. When considering the simulations in each context, calculate the effectiveness of correct classifications made by the algorithm in that context. This is done by comparing the NRMSE indexes of rejected and candidate solutions against the canonical analytical solution of the OCP, settling a certain error tolerance level to achieve a correctly classified solution as candidate. Solutions classified as candidates with indexes lower than the tolerance and the ones classified as rejected with indexes higher than the tolerance are considered to be correct classifications. The effectiveness in the context is given by η = NCorrect Classified Solutions NTotal of Solutions ; 4. The corroboration will occur when, after calculating the effectiveness of the algorithm in various contexts, the Null Hypothesis (H0) defined in Section 2.3 is corroborated by a statistical method, that is, the estimation of the median of the effectiveness is equal to or greater than 70%; and 2In fact, it is recommended to choose several OCPs, as done in this work. 38 Chapter 2. CONSIDERATIONS ABOUT THE METHODOLOGY 5. The rebuttal will occur when, after systematic and consistent testing, the Null Hypothesis (H0) gets refuted, raising the evidence of the Alternative Hypothesis (H1), with the median estimation of the effectiveness of the algorithm lower than 70%. It should be highlighted that the exclusion of solutions performed by the solutions classification method is a technical non-subjective procedure. On the other hand, defining whether a recommended solution is acceptable or not is a subjective evaluation procedure and may not be used in the verification process. However, once a solution is obtained, and after its classification as candidate or rejected, its confrontation against the exact solution in its context via NRMSE is a quantitative assessment perfectly compatible with the objective criteria to calculate the effectiveness of the algorithm. 2.5 STATISTICAL ANALYSIS AND TOOLS TO VERIFY THE HYPOTHESIS The Statistics Science is primarily applied in biological and human sciences (TRIOLA, 1999) in the study of processes involving random variables, which are usually non-numerical. However, its application in engineering is also common (LEVINE, 2000). In this Thesis, statistical methods were used with the help of specific tools to estimate the effectiveness of the proposed algorithm applicability. It is not scope of this study to address, in depth, statistical concepts. However, a brief consideration about the problem of the statistical evaluation of an algorithm is necessary. Therefore, the following assumptions were considered: • With the parametric variation introduced in steps 2 and 3, as well as the possible existing solutions in a nonconvex OCP, there will probably be a diversity of numerical solutions as high as the nonconvexity characteristics, the greater the parametric variations. This diversity of solutions and the limitations of the automatic classification algorithm suggest a 2.5. Statistical Analysis and Tools to Verify the Hypothesis 39 distribution of the effectiveness of correct classification of solutions for various simulations in several contexts. • It is considered that the algorithm has limited effectiveness at the minimum of 0% and maximum of 100%. It is also supposed an unknown distribution shape, as well as unknown mean and median. • It is necessary to be cautious and to have computational experience to efficiently apply an algorithm. A definition of an unfeasible OCP or inappropriate and unreasonable choice of parameters may lead to biased results that make it difficult to use the algorithm with efficiency. On the other hand, it is expected that, with the appropriate configuration and definition of the OCP, the algorithm is effective in the proper classification of numerous solutions obtained in an automated manner, saving considerable user time and effort in the individual classification of solutions. Given the assumptions above, and specially because one cannot affirm about the shape of the distribution, the Sign Test was chosen, as described in (TRIOLA, 1999, p.321) and (GIBBONS, 2003), to estimate the median of the effectiveness. This method was chosen for being "non-parametric", therefore, it does not require Normal distribution. The method was applied by own coding in Python language and using the routine "signtest()" from the MATLAB R© Statistics Toolbox. In this Thesis, 4 nonconvex OCPs were studied and a total of 972 simulations were optimized and analyzed in an automated way. It was possible, with the help of a second statistical method, to infer more information about the effectiveness distribution of the algorithm application in these OCPs. The second method is known as Kruskal-Wallis Test, (DANIEL, 1990). This test is also non-parametric and does not assume Normal distribution. However, the method has the following requirements (TRIOLA, 1999, p. 329): 40 Chapter 2. CONSIDERATIONS ABOUT THE METHODOLOGY • Minimum of 3 samples (there are 4 studied OCPs); • The test verifies the null hypothesis in which the samples derive from the same population or identical populations. In numerical terms it means to have the same median and the same distribution; • Each sample must have a minimum of five observations. In this Thesis there were 9 observations of effectiveness (one for each context in each OCP); and • The different samples cannot have very different variances. With the results of the simulations, the Kruskal-Wallis Test was applied with the aid of three different instruments: 1. By the function “scipy.stats.mstats.kruskalwallis()” of Python language; 2. By the function “kruskalwallis()” from MATLAB R© Statistics Toolbox; and 3. By using the application IBM R© SSPS R© Statistics, which results are presented in appendix C. The first two procedures for the Kruskal-Wallis Test resulted in the non-refuttal the Null Hypothesis (H0) in which the four samples had the same median. The latter presented as additional results an estimation of this median and the confirmation that the samples present the same distribution. The results, as detailed in Chapter 4, THE HYPOTHESIS VERIFICATION, were consistent. The next chapter presents the contribution of the Thesis, primarily by the proposition of an algorithm to treat nonconvex OCPs and, secondarily, by the structure ASAF - Automated Simulation and Analysis Framework, a set of Python scripts for automated simulations of OCPs. 41 3 CONTRIBUTIONS Here, the main contributions of this doctoral Thesis are presented. Initially, the details of the algorithm to treat nonconvex OCPs are presented as the main result of the Thesis. Next, there is a description of JModelica, IPOPT and other tools used as support for the research. Last, the ASAF - Automated Simulation and Analysis Framework is presented as secondary result. 3.1 ALGORITHM TO TREAT NONCONVEX OPTIMAL CONTROL PROBLEMS The treatment of nonconvex OCPs by automated numerical simulations can be carried out following the steps defined in Algorithm 1: Algorithm 1: Algorithm to Treat Nonconvex OCPs Data: 1)A problem that generates an OCP model; 2)Candidate solution requirements (α ,ε ,ρ)a; Result: 1)Classification of numerically viable solutions (feasible OCP); 2)Recommended solution (feasible OCP); repeat The Modeling Phase; 1 model the system, as well as the performance index, with fidelity to the real needs of the real problem and within acceptable domainsb: (DXi = [xi,L,xi,U ], i = 1...n, Du j = [u j,L,u j,U ], j = 1...m); ; 42 Chapter 3. CONTRIBUTIONS 2 make use of parameterization resources of the modeling and simulation tools to change the conditions of the model (system and performance index) within the limits defined in step 1; The Phase of Automated Obtaining of Solutions; 3 solve, in an automated manner, the OCP for parametric variation allowed by the modeling language, by the simulation and optimization tools and by the computational cost; The Phase of Automated Classification of Solutions; 4 analyze the results in an automated way for the previous parametric variation, define and distinguish regions, depending on the values of the parameters, in which the results of the OCP solutions represent a similar behavior; 5 exclude solutions, although they are numerically viable, that are inadequate from the practical applicability point of view, e.g., trivially null solutions; 6 exclude solutions that clearly do not comply with the restrictions of the classic theorems of existence (convexityc, non-decreasing monotonicityd and smooth growth of the performance index). The remainder solutions are called candidate solutions; foreach Solution Sν ∈ Solution Set S do r[0]← 0;mon← 0;convx← 0;convu← 0; for i← 0 to l−1 (where l = number of discretization points) do r[i+1] = dJ dt [i+1]/abs(J[i+1]) ; if (r[i+1]<0)||(r[i+1]>ρ) then mon++; for j← 1 to n do if (x j[i+1]− x j,L < ε)||(x j,U − x j[i+1]< ε) then convx ++; for j← 1 to m do if (u j[i+1]−u j,L < ε)||(u j,U −u j[i+1]< ε) then convu ++; if (mon∗100/l > α)||(convx ∗100/(n∗ l)> α)||(convu ∗100/(m∗ l)> α) then Rejected Solution; else Candidate Solution; The Choosing of the Recommended Solution; 7 among the remaining regions, choose the region that best suits the reality of the problem modeled by the OCP. The best candidate solution in this region, the recommended for this region, will be the recommended solution for the OCP; until The Remodeling Phase is exhausted with infeasible OCP for the User’s requirements or the recommended solution found is acceptable; aThese parameters are essential, however, they are not exhaustive in the implementation of the algorithm. Eventually, the user may reset them and adapt them to the characteristics of his or her definitions of simulations and contexts. bDo not to worry in advance if the resulting model does not comply with the conditions of linearity, convexity or monotonicity criteria of classical existence theorems. cThis is done by observing that the solution point (t,x) should be interior-point (CESARI, 1983, p. 107). In practical terms, the control and state variables cannot reach the boundaries of the OCP. dThis is done by direct inspection of performance index behavior. 3.1. Algorithm to Treat Nonconvex OCPs 43 In fact, the expression "interior-point"1 is related to a specific solution point in time. Perhaps for all solution points a better expression could be "Interior-Trajectory". Remark 3.1. The symbologies [•]L and [•]U stand for lower and upper limits, respectively. Remark 3.2. The symbology conv[•] denotes a nonconvexity counter. Remark 3.3. The variable mon denotes a non-compliant monotonicity counter. The notations in the algorithm, state vector x and control vector u, are used according the Optimal Control Theory and the definition of OCP presented in Appendix A. In Algorithm 1, the value of variable n represents the dimension of state vector x, the value of variable m represents the dimension of control vector u and the value of variable l represents the length of the numerical solutions vectors defined by the time discretization. Therefore, the variable l corresponds to the total number of collocation points. In this Thesis, for the practical implementation of the Algorithm 1, the following parameters were used: Table 3 – Parameters for the Solution Classification as Candidate Parameter Description Value α Maximum acceptable index for non-compliance 0,3(%) (nonconvexity, decreasing monotonicity, not slow growth of the performance index) ε Minimum acceptable slack to the limits of the variables 1e-4 (indicates the saturation of the variable) ρ Maximum acceptable growth rate of the Performance Index 0,1 Source: The author. The sequence of actions taken in the use of the algorithm is shown in Figure 2. 1The understanding of this expression here is similar, however, not the same, as in a classical static optimization, in a Nonlinear Program (NLP). It is similar to force the constraint that solution should not present active inequality constraints (LUENBERGER; YE, 2008, p. 322), more precisely, the solution should not be at the frontier of the feasibility region (LUENBERGER; YE, 2008, p. 322) 44 Chapter 3. CONTRIBUTIONS Fi gu re 2 – Se qu en ce fo rt he A pp lic at io n of th e A lg or ith m .S ou rc e: T he au th or . 3.1. Algorithm to Treat Nonconvex OCPs 45 It is understood that numerical simulations algorithms, in which parameters present some randomness or arbitrary choice2, carry this sensitivity to the possible solutions that were reached. Therefore, depending on the definition of the problem, the choices of the parameters, and the methods used in numerical simulation, it may not converge, converge to a specific solution and, once the convergence is reached, the solution may be acceptable or not. The proposed algorithm presents a quantitative dimension during the phase when the candidate solutions are defined. However, the choice of acceptable solution among the various recommended ones is made in a subjective manner by the user. The quantitative assessment of the solutions can be easily incorporated in the automation process, however, the qualitative phase, consisted by the definition about the acceptability or non-acceptability of the recommended solution, cannot and should not be automated, since it is a result of the judgment of qualitative merit made by the human element who made the initial parametric choices. It was verified that the initial freedom the user has to make parametric choices can affect three very specific and delimited groups of parameters: 2As in the mutation factor in genetic algorithms, or the definition of the network structure, number of neurons, layers and feedback of an artificial neural network. 46 Chapter 3. CONTRIBUTIONS Definition 3.1. Type-I Parameters Group: It constitutes the group of parameters that affect the equations that model the system in study, including those related to the law of control3. Definition 3.2. Type-II Parameters Group: They are related to the essential and specific parameters of simulation and optimization tools4, such as number of discretization elements, steps, tolerances and stopping criteria. The definition of algorithms and heuristics used by the tools can be found in this group. Definition 3.3. Type-III Parameters Group: Contains the parameters determined by the characteristics of the OCP Design and Constraints. They are restrictions of lower or upper limits to the state variables, to the control variables, the initial conditions of these variables, if fixed or free, the parameters that affect the performance index and, eventually, other, if they exist5,6. Starting time and Ending time of the simulations are examples of parameters of this group. It should be noted that the parameters variation procedure performed in the initial phase of the implementation of the algorithm shown in Figure 2, allows a diversity of solutions and, in most situations, various contexts. In this Thesis and specifically to perform the verification of the hypothesis, context is defined as: Definition 3.4. Context: 3Without loss of generality, but specifically to the tools chosen in this Thesis, the model parameters of the system are the ones that affect the section of equations defined after the keyword equation of Modelica language, at the file of the system model. 4Actually, all the parameters, including the type-I and type-III, are set in the tools, but one can understand as type-II parameters the ones that do not affect the system model neither the design of the OCP, they are specific and unique of the tools and their algorithms, such as tolerances, steps and stopping criteria. 5Eventually the control may require that a variable have limits for its derivative. 6Similarly, for the tools used based on modeling in Modelica language, typical examples are the constraints inserted in the model after the keyword constraint. 3.1. Algorithm to Treat Nonconvex OCPs 47 It is a specific and unique choice of values for type-I and type-III parameters. All parameters that are able to modify the general settings of the OCP, either by the change in the system model, via state equations with the law of control, or by change in the objective function, or by restriction or delimitation of the possible solution sets of state and control variables. The definition of Context is required to make the objective comparison of solutions consistent, regardless if they are rejected or candidate, under a same numeric scale, such as the NRMSE. It is a definition necessary only in the process of validation of the hypothesis. In the general scope of the algorithm application, this definition is unnecessary because there is no formal application of metric scale to candidate solutions. Figure 3 presents a graphical visualization of the use of the algorithm and the delimitation of contexts and regions. It is appropriate to note that the proposed algorithm has a generic characteristic and does not depend on a specifically determined tool or context for its application. In fact, the final choice of the recommended solution is a qualitative judgment of merit by the user, in consonance with the parametric choices, that mirror his or her real requirements for the solution of the OCP. However, in the scope of verification of the Hypothesis, the definition of the concept of Context is necessary and essential. Another important note is that the term region used in the algorithm is more generic and is related to the perspective the user has about the results of the simulations. In fact, a region could be a group of solutions, a unique and isolated context or even a group of contexts. Next, the details about the integrated simulation environment that was used to perform the optimizations that supported the Thesis are presented. 48 Chapter 3. CONTRIBUTIONS Figure 3 – Description of the Process of Algorithm Application. Source: The author. 3.2 THE INTEGRATED AND AUTOMATED ENVIRONMENT FOR MODELING, SIMULATION AND OPTIMIZATION The proposed algorithm tries to make better use of the parameterization capabilities and automation of modern computational tools. It is well known that nonlinear and nonconvex systems, in most cases, have different solutions and behaviors due to parameterization and initial conditions. To analyze such systems, researchers use classical instruments such as Attractors, Basins of Attraction and Bifurcation, Phase Plane Analysis and Phase Portraits, (THOMPSON; STEWART, 2002). However, (SLOTINE; LI, 1991, p.17) affirms the natural complexity of the graphic study of systems of order 3 and greater. Thus, the application of the proposed algorithm demanded the implementation of a computational environment that deals with this matter by holding hundreds, even thousands, 3.2. The Integrated and Automated Environment for Modeling, Simulation and Optimization 49 of optimizations and simulations with the generation of huge amount of raw data for further analysis. Then, it was carried out the development of a code in Python scripts specific to the Thesis, defined by the ASAF-Automated Simulation and Analysis Framework structure, presented in Section 3.3. The ASAF-Automated Simulation and Analysis Framework has its execution based on the following tools. 3.2.1 JModelica.org The JModelica.org tool (http://jmodelica.org/) (ÅKESSON, 2010) is an integrated environment of modeling, simulation and optimization of dynamic systems. This is the result of a research held by Department of Automatic Control at the University of Lund, Sweden, that currently has its development and maintenance made by the Swedish company Modelon AB. A description of the JModelica.org integration with other tools, especially the flow of data and commands, is shown in Figure 4. This simulation environment was deployed following the tutorials provided by maintainers of the tools, (ÅKESSON, 2010), (ANDERSSON, 2011; MAGNUSSON; ÅKESSON, 2012), (WÄCHTER; BIEGLER, 2006), in a more evolved form than what was presented in Sections 4.2 to 4.5 of (RESENDE, 2015). The main update was regarding the operating system, which currently is the Ubuntu LTS Desktop 16.04, corresponding to the version with support available until the year 2021 (Long-Term Support). The current versions of the tools are presented in the Table 4. Through the observation of Figure 4, one can verify that the user calls directly the JModelica.org routines, which communicates with other tools through internal routines. For example, the user does not interact directly with the IPOPT optimizer, but one can, by calling 50 Chapter 3. CONTRIBUTIONS Figure 4 – Automated Structure to Solve Nonlinear and Nonconvex OCPs. Source: The author. the JModelica.org routine, set the tolerances and other parameters of the optimizer. Analogously, after the optimization, IPOPT does not report the results directly to the user by saving the values of variables, but rather through calls to routines and objects of the JModelica.org. The main feature of JModelica.org in this environment is to interpret an OCP coded in Modelica language, translating the OCP to a Nonlinear Program (NLP) to be optimized by IPOPT. 3.2.2 IPOPT The IPOPT tool, (WÄCHTER; BIEGLER, 2006), is compiled separately from the JModelica.org, however, in this environment, its operation is integrated to JModelica.org. Once these tools are compiled, the user only has to set the optimization options, in case of different 3.2. The Integrated and Automated Environment for Modeling, Simulation and Optimization 51 values from the default ones, and execute the optimization. The interesting part about this integrated environment is that all options available in the IPOPT User Manual can be set via JModelica.org: The ones that affect the stopping mechanisms, the speed and convergence of iterations, including options concerned to the barrier method, KKT conditions and others. It is evidenced that the good choice of these options can determine the success7 and the speed of the optimization. It is important to emphasize that the IPOPT performs a purely static optimization. It is up to the JModelica.org to transform the intrinsic dynamism of the original OCP to a static NLP, as described in appendix A. The IPOPT has an integrated library called sIPOPT that allows the execution of sensitivity analysis. 3.2.3 HSL Routines The HSL routines, according to (Science and Technology Facilities Council, 2015), constitute a set of solvers for large-scale linear problems. These routines are provided under the “HSL ACADEMIC LICENSE VERSION 1.1 JANUARY 2011” license, and can be compiled in an integrated manner to IPOPT or dynamically loaded. Since there are several routines, once compiled, the usage is made by setting the option “solver” of the IPOPT. 3.2.4 Python and other tools The JModelica.org, IPOPT and HSL routines constitute the core of the modeling, optimization and simulation environment. However, other tools are required both in the installation process of the environment as well as in its operation. The Python language interpreter itself must be installed, (ROSSUM, 1995). In addition, 7An important factor that directly affects the optimization is the consistent definition of the OCP constraints. Eventually, the IPOPT will not converge due to the simple fact that the OCP has been formulated so restrictively that it became an infeasible problem. 52 Chapter 3. CONTRIBUTIONS it is necessary to have a compiler for the C programming language and the Java Development Kit (JDK) installed for the compilation and operation of the JModelica.org, IPOPT and HSL tools. If the environment is Microso f t R©Windows R©, one can install the binaries downloaded from the JModelica.org Web site. This Thesis, however, has used the Linux environment8 with the tools in Table 4. Table 4 – Tools Used in this Thesis. Tool Version License JModelica.org 1.17, 2.0, 2.1 and 2.2 GPL v3 IPOPT 3.12 Eclipse Public License - v 1.0 HSL 2014.01.10 HSL ACADEMIC LICENSE VERSION 1.1 JANUARY 2011 Python 2.7.12 PYTHON SOFTWARE FOUNDATION LICENSE VERSION 2 Java (OpenJDK) 1.8.0_131 GPL v2+ gcc 5.4.0 GPL v3 Source: Data provided publicly in the tools documentation. 3.3 ASAF - AUTOMATED SIMULATION AND ANALYSIS FRAMEWORK Here, the specific Python and Modelica codes developed for the Thesis are presented as secondary result of the research. The ASAF - Automated Simulation and Analysis Framework is limited to the Python scripts and the Modelica model required to perform the optimizations and simulations. The scripts that comprise the ASAF structure are: • Pack.mop - This is an entry file to the JModelica.oirg, containing the system modeled in Modelica language; • asaf_config_.py - This is the configuration file that sets the main parameters for simulations and optimizations; • asaf_pgm.py- This is the script that is responsible to perform the simulations and optimizations; 8The operating system used in this Thesis was the Ubuntu LTS Desktop 16.04 3.3. ASAF - Automated Simulation and Analysis Framework 53 • asaf.py - This is the file with the code of Python classes that are specific to the ASAF structure. It is a generic code used to treat dynamic systems; • asaf_global.py - This file contains global functions called by the main script; and • asaf_graph.py - This file contains graphics functions necessary to get the results of the analysis, as plots of state and control signals. All these files constitute a framework for studying nonlinear dynamical systems based on the environment described in Section 3.2. To use this structure, users have to edit only two files: Pack.mop and asaf_config_.py. The procedure to use this environment to study a specific system starts by modeling the system in Modelica language directly into the file Pack.mop. Next, comes the editing of the asaf_config_.py configuration file according to the model, to set the parameters for simulations/optimizations. The script asaf_config_.py is executed and all the configuration variables are stored in the file asaf.dict, which will be read by the main script asaf_pgm.py. After that, the main script asaf_pgm.py is executed. It will generate all the data and figures resulting from the simulations/optimizations. The numerical results are written in files type “.h5” (Hierarchical Data Format - HDF) and the figures are stored in graphic formats defined in the configuration (e.g. ".jpg" and ".eps"). An important consequence of the use of integrated and automated environment is that, once installed the tools presented in Table 4, the operation of the ASAF structure, with the use of the respective configuration and the Modelica model files for a given system/OCP, allows total and complete reproduction of the numerical results and graphics, not only of this Thesis, but for any OCP for which one has stored the model and configuration files. 54 Chapter 3. CONTRIBUTIONS 3.3.1 Main Features of the ASAF Structure The main features of the ASAF application can be used by setting two of its most important configuration items: The Running_Type and the Running_Procedure variables. The current options for Running_Type are: • Optimization: to perform optimizations - in this case the system model, jointly defined by model_name and Opt_Label variables, must be coded with Optimica extension of Modelica language; • Simulation: to perform simulations; • Initialization: to attempt to find equilibrium points of a system; and • Analysis: to perform analysis, generating graphics and animations from raw data calculated in previous simulations or optimizations. The current options for Running_Procedure are: • Default: to carry out an analysis, simulation or optimization of a model without parametric variation; • Cases: to carry out an analysis, optimizations or simulations of a model with parametric variation corresponding to specific cases of vectors of different parameters defined in a group; • Parameters: to carry out an analysis, optimizations or simulations of a model with parametric variation corresponding to specific cases of vectors of predefined parameters. Each parameter has its vector of values and every simulation or optimization occurs after a change, from the model with default values, of a single parameter value at a time, therefore covering all the vector of values of each parameter; • Phase_Portrait: for the generation of Phase Portrait with simulations of a model; 3.3. ASAF - Automated Simulation and Analysis Framework 55 • Mixed_Cases_and_Parameters: jointly Cases and Parameters; • Mixed_Phase_Portrait_and_Parameters: jointly Phase_Portrait and Parameters; and • Mixed_Phase_Portrait_and_Cases: jointly Phase_Portrait and Cases. The features described above are mature, stable and can be fully applied in the studies of dynamical systems and solutions of OCPs. However, other features can and will be under development, such as Basins of Attraction, Poincaré maps, Lyapunov Exponents calculation and 0-1 Test (GOTTWALD; MELBOURNE, 2009). A project of graphical user interface for the ASAF structure is under development: 56 Chapter 3. CONTRIBUTIONS Fi gu re 5 – D ev el op m en to ft he G ra ph ic al U se rI nt er fa ce .S ou rc e: T he au th or 3.3. ASAF - Automated Simulation and Analysis Framework 57 Figure 5 shows the configuration used for simulations with OCP#4. The development of this easiness to the user aims to minimize the chances of problems and mistakes in parametric definitions, allowing direct and objective choice of parameters, as well as to ensure the proper registration on file for all Simulations/Optimizations performed. The importance of this feature is evidenced by the fact that the parametric definition (such as the steps and barrier parameter) is essential for an algorithm to achieve numerical stability (PARKER; CHUA, 1991, p. 97) and can determine whether or not the convergence of iterations will occur, and if the solution will be achieved. As application of this structure in the study of dynamical systems, in (RESENDE, 2017) it was presented numerical simulations of a generalized Lorenz system with four dimensions and hyperchaotic behavior, according to the model presented in (MACEK; STRUMIK, 2010) and (MACEK; STRUMIK, 2014). Next, in the following chapter it will be presented the entire procedure performed for an initial verification of the hypothesis. 59 4 THE HYPOTHESIS VERIFICATION The main caution that should be taken in an attempt to test the Null Hypothesis (H0) using the method defined in Section 2.4, is related to the choice of OCPs to be used. The chosen OCPs must contain two very specific features: It must be “non-compliant” with the classical theory, which means that it must be nonconvex, and it must have an exact solution known. These conditions for the OCPs that are candidates to be used in the verification of the hypothesis makes it very difficult to find appropriate OCPs for the examination. The following sections present the OCPs chosen for the Hypothesis Verification. The analytical solutions, necessary to establish a comparison reference for the numerical solutions provided by the application of the algorithm, are also presented. 4.1 FIRST OCP USED IN THE HYPOTHESIS VERIFICATION After initial research, the following problem was found in (MEZIAT, 2007, p.156): minimize: J [xxx,uuu, t] = ∫ 1 0 (x−10∗ t)2 dt (1) subject to: ẋ = u2−u∗ x+ x (2) x(0) = 0 (3) 60 Chapter 4. THE HYPOTHESIS VERIFICATION By observing the OCP defined above, one may verify that it is nonlinear and nonconvex. It can be verified by direct inspection of the performance index, that it is possible to infer the state optimal trajectory as xxx∗ ≡ x(t) = 10∗ t. By substituting the expression for x(t) in the equation (2), the following possibilities for the control vector uuu∗ are obtained: u∗1(t) = 10∗ t− √ 100∗ t2−4∗ (10∗ t−10) 2 (4) u∗2(t) = 10∗ t + √ 100∗ t2−4∗ (10∗ t−10) 2 (5) The curves of the respective solutions are presented, respectively in Figures 6 and 7: 0.0 0.2 0.4 0.6 0.8 1.0 t 3 4 5 6 7 8 9 10 u ∗ 1 First Analytical Control Solution Figure 6 – First Analytical Solution of OCP#1. Source: The author. 4.1. First OCP Used in the Hypothesis Verification 61 0.0 0.2 0.4 0.6 0.8 1.0 t −3.5 −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 u ∗ 2 Second Analytical Control Solution Figure 7 – Second Analytical Solution of OCP#1. Source: The author. 62 Chapter 4. THE HYPOTHESIS VERIFICATION 4.2 SECOND OCP USED IN THE HYPOTHESIS VERIFICATION After research, the following problem was found in (MEZIAT, 2007, p.157): minimize: J [xxx,uuu, t] = ∫ 10 0 ( x− t2)2 dt (6) subject to: ẋ = u2−u∗ x+ x (7) x(0) = 0 (8) Again, it is possible to infer the state optimal trajectory as xxx∗≡ x(t)= t2 and the following possibilities for the control vector uuu∗: u∗1(t) = t2 + √ t4−4∗ (t2−2∗ t) 2 (9) u∗2(t) = t2− √ t4−4∗ (t2−2∗ t) 2 (10) The curves of the solutions are: 4.2. Second OCP Used in the Hypothesis Verification 63 0 2 4 6 8 10 t −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 u ∗ 1 First Analytical Control Solution Figure 8 – First Analytical Solution for OCP#2. Source: The author. 0 2 4 6 8 10 t 0 20 40 60 80 100 u ∗ 2 Second Analytical Control Solution Figure 9 – Second Analytical Solution for OCP#2. Source: The author. 64 Chapter 4. THE HYPOTHESIS VERIFICATION 4.3 THIRD OCP USED IN THE HYPOTHESIS VERIFICATION With the aim to amplify the possibilities of verifying the hypothesis, a third OCP was chosen from (MEZIAT, 2007, p. 158): minimize: J [xxx,uuu, t] = ∫ 1 0 ( (x− t)2 +(y− t)2)dt (11) subject to: ẋ = u (12) ẏ = (1−u2)2 + x2 (13) x(0) = 0; y(0) = 0 (14) Once again, it can be observed by direct inspection of the performance index that it is possible to infer the optimal trajectory as xxx∗ ≡ (x(t),y(t)) = (t, t). This does not mean that this ideal trajectory is feasible by the system. For this OCP, however, even if it is possible to find the state ideal solution, an analytical solution for the control vector uuu∗ was not obtained. 4.4 FOURTH OCP USED IN THE HYPOTHESIS VERIFICATION Finally, the problem presented in (MEZIAT, 2007, p. 160) was chosen as fourth OCP: minimize: J [xxx,uuu, t] = ∫ 1 0 ( (x− t)2 +(y− t)2)dt (15) subject to: ẋ = u (16) ẏ = (1−u)2 ∗ (2−u)2 + x2 (17) x(0) = 0; y(0) = 0 (18) 4.5. Results 65 Once again, it is possible to infer the optimal trajectory as xxx∗ ≡ (x(t),y(t)) = (t, t). This does not mean that this ideal trajectory is feasible by the system. For this OCP, as in the previous OCP, an analytical solution for the control vector uuu∗ was not obtained. 4.5 RESULTS Table 5 presents a summary of the canonical solutions for the used OCPs. Table 5 – Summary of the Canonical Solutions of the Used OCPs OCP #1 OCP #2 OCP #3 OCP #4 u∗1 10∗t− √ 100∗t2−4∗(10∗t−10) 2 t2+ √ t4−4∗(t2−2∗t) 2 - - u∗2 10∗t+ √ 100∗t2−4∗(10∗t−10) 2 t2− √ t4−4∗(t2−2∗t) 2 - - x∗ 10∗ t t2 t t y∗ - - t t Source: The author. As an attempt to verify the research hypothesis, optimizations with the default parameters of the IPOPT optimizer defined in version 3.12 were performed, except for the changes defined in Table 6: Table 6 – Optimization Parameters Used in All the Optimizations Parameter Tool Description Value n_c JModelica.org As defined in List of Symbols 10 linear_solver IPOPTa Linear Solver ma57 max_iterb IPOPT Maximum number of iterations 3e4 mu_maxb IPOPT Maximum value for the barrier parameter 2 mu_minb IPOPT Minimum value for the barrier parameter 1e-300 a According to Section C, Options Reference, from the User’s Guide Introduction to Ipopt: A tutorial for downloading, installing, and using Ipopt,of the IPOPT 3.12 tool. b Has effect only if the option “mu_strategy” chosen is “adaptive”. Source: The author. The definition of each context was made by choosing two parameters, setting the duo: C = {PC1,PC2} (19) 66 Chapter 4. THE HYPOTHESIS VERIFICATION Where PC1 and PC2 were determined by the conditions presented in Table 7. Table 7 – Parameters Used in the Definitions of the Contexts OCPa Index Parameter (s) to Set Possible Values #1 #2 #3 #4 I [2, 12] [-2, 2] [0, 2] [-2, 10] PC1 Interval [umin,umax] II [-5, 0] [-1, 110] [0,7, 5] [-5, 1] III [-2, 8] [0, 64] [-1, 0,7] [1, 10] I c 0 0 0 0 PC2 periodb II 0,05 0,5 0,05 0.05 III 0,1 1 0,1 0,1 a Specific values for each OCP b A parameter that sets the “blocking_factors” option of the JModelica.org tool, used for piecewise constant controls. c With “period” = 0, here is no restriction of constant controls in periods. Source: The author Therefore, with the possibilities of adjustments on Table 7, nine contexts are defined for each OCP. The definition of each simulation was made by choosing four parameters, setting the tuple: S = {PS1,PS2,PS3,PS4} (20) Where each PSi was determined by Table 8. Therefore, with the possibilities of adjustments on Table 8, twenty-seven different simulations are defined, that remain the same for each context of each OCP. Table 9 presents a summary of the scope of the simulations. The parameters of Table 6 were kept the same in all the optimizations. The parameters in Tables 7 and 8 were adjusted accordingly for each specific simulation in a given context. Next, the obtained results are presented. 4.5. Optimizations Results 67 Table 8 – Parameters Used in the Definitions of the Simulations Index Parameter to Set Possible Values Set I 40 PS1 n_ea II 200 III 1000 I adaptive PS2 mu_strategyb II monotone I kkt-error PS3 adaptive_mu_globalizationc II obj-constr-filter I safer-min-dual-infeas PS4 alpha_for_yd II primal III full a The parameter that sets the option “n_e” of the JModelica.org tool, used to determine the number of elements in the discretization. b The parameter that sets the option “mu_strategy” of the IPOPT tool. c The parameter that sets the option “adaptive_mu_globalization” of IPOPT tool. Has effect only if the option “mu_strategy” chosen is “adaptive”. d The parameter that sets the option “alpha_for_y”of the IPOPT tool. Source: The author. Table 9 – Scope of Simulations Total of OCPs Contexts per OCP Simulations per Context Total of Simulations 4 9 27 972 Source: The author. 4.5.1 Optimizations Results for OCP#1 Here, the results of the optimizations are presented: Table 10 – Optimizations Results for OCP#1 Context Candidates Rejected Not Solved Correct Classifications Effectiveness (%) A 21 6 0 21 78 B 25 2 0 27 100 C 27 0 0 27 100 J 27 0 0 18 67 K 27 0 0 27 100 L 27 0 0 0 0 S 13 9 5 15 68 T 9 16 2 17 68 U 14 13 0 25 93 Average - - - - 74.89 Source: The author. 68 Chapter 4. THE HYPOTHESIS VERIFICATION 4.5.1.1 Analysis of the Optimizations Results for OCP#1 It is observed, in context A, that because of the definition of the control variable limits, the results were similar to the first canonical solution as shown in the figures of the simulation A-1: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 3 4 5 6 7 8 9 10 u Control - NRMSE=0.20386% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 10 – Candidate Solution for OCP#1 - Context A: Control. Source: The author. 4.5. Optimizations Results 69 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.19835% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 11 – Candidate Solution for OCP#1 - Context A: State. Source: The author. 70 Chapter 4. THE HYPOTHESIS VERIFICATION Diversely, the solution was rejected in the simulation A-4. The rejection did not occur due to an inadequacy of the found state, but rather due to the high NRMSE index of the control variable: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 2 4 6 8 10 u Control - NRMSE=2.83217% (Lower Sat= 1.50%; Upper Sat= 0.00%; Total Sat= 1.50%) u u1_ref u2_ref Figure 12 – Rejected Solution for OCP#1 - Context A: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.19849% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 13 – Rejected Solution for OCP#1 - Context A: State. Source: The author. 4.5. Optimizations Results 71 The solution A-17 was numerically classified as a candidate by the algorithm. However, it is observed in Figure 14 that the solution is not qualitatively interesting. This wrong classification was correctly identified in the verification of the hypothesis, when performing the automated confrontation against the ideal analytical solution. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 3 4 5 6 7 8 9 10 u Control - NRMSE=4.85196% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 14 – Solution Improperly Classified as Candidate for OCP#1 - Context A: Control. Source: The author. 72 Chapter 4. THE HYPOTHESIS VERIFICATION 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.03925% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 15 – Solution Improperly Classified as Candidate for OCP#1- Context A: State. Source: The author. 4.5. Optimizations Results 73 Similar situations have occurred in all other simulations of this context. Context B also presented solutions that leaned towards to the first canonical solution, respecting the piecewise control restriction, as shown by the figures of simulation B-2, as follows: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 3 4 5 6 7 8 9 10 u Control - NRMSE=1.48820% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 16 – Candidate Solution for OCP#1 - Context B: Control. Source: The author. 74 Chapter 4. THE HYPOTHESIS VERIFICATION 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.06508% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 17 – Candidate Solution for OCP#1 - Context B: State. Source: The author. 4.5. Optimizations Results 75 Simulation B-1 was rejected due to the high level of saturation index of the control variable: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 2 4 6 8 10 u Control - NRMSE=6.07036% (Lower Sat= 5.00%; Upper Sat= 0.00%; Total Sat= 5.00%) u u1_ref u2_ref Figure 18 – Rejected Solution for OCP#1 - Context B: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.73758% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 19 – Rejected Solution for OCP#1 - Context B: State. Source: The author. 76 Chapter 4. THE HYPOTHESIS VERIFICATION Context C also presented solutions that leaned towards to the first canonical solution: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 3 4 5 6 7 8 9 10 u Control - NRMSE=2.97513% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 20 – Candidate Solution for OCP#1 - Context C: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.29105% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 21 – Candidate Solution for OCP#1 - Context C: State. Source: The author. 4.5. Optimizations Results 77 Without presenting further figures, a similar phenomenon was observed for the J, K and L contexts, this time, in the range of the second canonical solution. However, something really interesting happened in contexts S, T and U, where the range of the control partially encompassed both canonical solutions. Thus, it was found in many simulations the occurrence of numerical solution with oscillation points between the canonical solutions, as shown in the following figures. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] −2 0 2 4 6 u C ntr l - NRMSE=0.08384% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 22 – Candidate Solution for OCP#1 - Context S: Control. Source: The author. 78 Chapter 4. THE HYPOTHESIS VERIFICATION 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.19835% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 23 – Candidate Solution for OCP#1 - Context S: State. Source: The author. 4.5. Optimizations Results 79 0.0 0.2 0.4 0.6 0.8 1.0 t[s] −2 −1 0 1 2 3 4 5 u Contro - NRMSE=3.57651% (Lower Sat= 6.75%; Upper Sat= 0.00%; Total Sat= 6.75%) u u1_ref u2_ref Figure 24 – Rejected Solution for OCP#1 - Context S: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.19777% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 25 – Rejected Solution for OCP#1 - Context S: Control. Source: The author. 80 Chapter 4. THE HYPOTHESIS VERIFICATION Similarly, it was observed in context T: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] −2 −1 0 1 2 3 4 u Contro - NRMSE=0.76849% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 26 – Candidate Solution for OCP#1 - Context T: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.19963% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 27 – Candidate Solution for OCP#1 - Context T: State. Source: The author. 4.5. Optimizations Results 81 0.0 0.2 0.4 0.6 0.8 1.0 t[s] −2 0 2 4 6 8 u C ntr l - NRMSE=6.08582% (Lower Sat= 0.00%; Upper Sat=25.00%; Total Sat=25.00%) u u1_ref u2_ref Figure 28 – Rejected Solution for OCP#1 - Context T: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=1.65046% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 29 – Rejected Solution for OCP#1 - Context T: State. Source: The author. 82 Chapter 4. THE HYPOTHESIS VERIFICATION Similarly, it was observed in context U: 0.0 0.2 0.4 0.6 0.8 1.0 t[s] −1 0 1 2 3 4 u Contro - NRMSE=1.59804% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) u u1_ref u2_ref Figure 30 – Candidate Solution for OCP#1 - Context U: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.22480% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 31 – Candidate Solution for OCP#1 - Context U: State. Source: The author. 4.5. Optimizations Results 83 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 u Control - NRMSE=3.78951% (Lower Sat= 0.00%; Upper Sat=20.00%; Total Sat=20.00%) u u1_ref u2_ref Figure 32 – Rejected Solution for OCP#1 - Context U: Control. Source: The author. 0.0 0.2 0.4 0.6 0.8 1.0 t[s] 0 2 4 6 8 10 x State Var 1 - NRMSE=0.94338% (Lower Sat= 0.00%; Upper Sat= 0.00%; Total Sat= 0.00%) x x_ref Figure 33 – Rejected Solution for OCP#1 - Context U: State. Source: The author. 84 Chapter 4. THE HYPOTHESIS VERIFICATION The simulations of OCPs from 2 to 4 presented a similar outcome. The detailed numerical report is presented in Appendix B. A brief summary of the results is presented in Table 11. Table 11 – Estimation of the Effectiveness of the Algorithm: Summary of the Results OCP# Effectiveness Average (%) JM1.17, JM2.0, JM2.1 and JM2.2 1 74.89 2 72.78 3 95.56 4 87.22 Global Average of the Effectiveness 82.61 Source: The author. The simulations were performed in four different versions of JModelica.org. The results were numerically the same in all versions. 4.6 THE HYPOTHESIS VERIFICATION BY STATISTICAL TESTS The Signal Test was performed in MATLAB R© with the following code: 4.6.1 The MATLAB R© Code for the Hypothesis Verification 1 %−−−−−−−−−−−−−−−−−−−−−−−−− 2 c l e a r a l l ; c l o s e a l l ; c l c ; 3 4 c l e a r a ; 5 a =[78 , 100 , 100 , 67 , 100 , 0 , 68 , 68 , 9 3 ; . . . 6 100 , 63 , 45 , 100 , 0 , 64 , 100 , 83 , 1 0 0 ; . . . 7 100 , 100 , 93 , 100 , 89 , 78 , 100 , 100 , 1 0 0 ; . . . 8 100 , 100 , 63 , 67 , 88 , 92 , 75 , 100 , 100 9 ] ; 10 11 a= round ( a ) 12 [ p , t b l , s t a t s ] = k r u s k a l w a l l i s ( a ’ ) 13 f i g u r e 4.6. The Hypothesis Verification by Statistical Tests 85 14 multcompare ( s t a t s ) 15 [ p1 , h1 , s t a t s ] = s i g n t e s t ( a ( 1 , : ) , 7 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ a l p h a ’ , 0 . 0 0 1 , ’ method ’ , ’ e x a c t ’ ) 16 [ p2 , h2 , s t a t s ] = s i g n t e s t ( a ( 2 , : ) , 7 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ a l p h a ’ , 0 . 0 0 1 , ’ method ’ , ’ e x a c t ’ ) 17 [ p3 , h3 , s t a t s ] = s i g n t e s t ( a ( 3 , : ) , 7 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ a l p h a ’ , 0 . 0 0 1 , ’ method ’ , ’ e x a c t ’ ) 18 [ p4 , h4 , s t a t s ] = s i g n t e s t ( a ( 4 , : ) , 7 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ a l p h a ’ , 0 . 0 0 1 , ’ method ’ , ’ e x a c t ’ ) 19 20 21 22 [ p1 , h1 , s t a t s ] = s i g n t e s t ( a ( 1 , : ) , 9 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ method ’ , ’ e x a c t ’ ) 23 [ p2 , h2 , s t a t s ] = s i g n t e s t ( a ( 2 , : ) , 9 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ method ’ , ’ e x a c t ’ ) 24 [ p3 , h3 , s t a t s ] = s i g n t e s t ( a ( 3 , : ) , 9 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ method ’ , ’ e x a c t ’ ) 25 [ p4 , h4 , s t a t s ] = s i g n t e s t ( a ( 4 , : ) , 9 0 , ’ t a i l ’ , ’ r i g h t ’ , ’ method ’ , ’ e x a c t ’ ) The Signal Test results corroborated the Null Hypothesis in which the median estimate of the effectiveness of the algorithm is equal to or greater than 70%. On the other hand, the Kruskal-Wallis Test, (DANIEL, 1990), corroborated the hypothesis that the four samples related to the four studied OCPs present the same median of effectiveness of the algorithm implementation, also being originated from the same population. Further detailing of the Kruskal-Wallis Test was achieved with the application IBM R© SPSS R© Statistics. The results can be found in Appendix C, where the median estimate for all samples was calculated at 93%, with significance of 0.05 (95%). Therefore, within the scope defined and based on the statistical tests performed, the hypothesis of this doctoral research is corroborated. The next chapter presents a brief discussion on the results of the research and its implications for the study and solution of nonconvex OCPs by means of automated numerical simulations. 87 5 DISCUSSION OF THE RESULTS The discussion of the results of this Thesis starts by evoking the experience lived by Umberto Eco, described in (ECO, 2007, p. 112)1, which resulted in his guidance on Scientific Humility: “the fact is that we must respectfully listen everyone, without ceasing to express value judgments...” In this Thesis, parallelly to the concept originally defined for bibliographic references, an attempt was made to apply Scientific Humility in the scope of the possible context contributions and their numerical solutions provided by the algorithm. Hence, one should not exclude in advance the project of possible contexts and, consequently, their numerical solutions for a given OCP. It is recommended that, for the efficient application of the proposed algorithm, one should vary parameters, define various contexts without bias, seek solutions in an automated manner, classify solutions with criteria, reject the inappropriate and, among the remainder, choose a recommended solution, if acceptable. With the continuous advancement of computational capacity and metaprogramming techniques, the approach by automated numerical simulations is an important resource that can, at first, perform numerous parametric variations to generate a diversity of numerical solutions for non-compliant problems, that cannot be solved by conventional methods. 1Text in Portuguese 88 Chapter 5. DISCUSSION OF THE RESULTS This Thesis explores exactly these new computational possibilities and features a criteria and a rational method, without excluding a context or a solution in advance, to provide a recommended solution to a nonconvex OCP, object of the algorithm. The proposed algorithm is generic enough to not depend on current computational tools and can be applied with other tools and computational techniques. Parallel computation or analysis of dynamical systems with the aid of virtual reality to study the results by 3D visualization, could be exploited jointly to the use of the algorithm. In the future, emerging technologies such as quantum computing may become viable and open new frontiers for the application of this algorithm. In fact, (BERNSTEIN, 2009, p. 19) states that the interest in quantum computation is least related to the algorithms and programs, but rather mostly to the ability of processing of data in a parallel way. Thus, the new technologies, whether they are conventional or emerging, point to a increasingly lower computational cost to perform automated numerical simulations and the respective processing of data. The algorithm of Section 3.1 stands out with the main result of the Thesis, however, the computational structure presented in Section 3.3 represents a very important secondary accessory contribution, since it made it possible to perform optimizations in an automated manner that, for now, corroborates the hypothesis. Usually, the academic study of dynamical systems encompasses the specific coding to solve each problem that is presented. The framework developed in this Thesis has the functionality to treat OCPs in a generic way and with very little coding work. The tools used by the framework have free or academic licenses and the user has full access to all features of the programs. The framework and tools used are updated to the most recent versions released by the developers. Any researcher who install the tools following their tutorials, can achieve the complete reproduction of all results of the simulations and optimizations by the automated execution of 89 the configuration codes specifically developed for this Thesis. In this Thesis, it is observed that in every simulation process, classifications of solutions and even of hypothesis verification, there was a high degree of automation and parameterization. It was considered an error threshold allowed for the NRMSE for a candidate solution over the ideal solution. In contexts without constant controls restrictions (PC1 = I), this threshold was 1%, in intermediate-constrained contexts (PC1 = II), this threshold was 2%, and in other contexts (PC1 = III), this threshold was 3%. Parametric variations which, at the stage of OCP modeling and optimizations are the origin of the diversity of the numerical solutions found, at the stage of classifications and verification of the hypothesis can lead to biases that may interfere in the results. Thus, it is important to have criteria to define them, especially the analysis parameters. It is also noted that the strictly numerical automated evaluation can differ from a qualitative individual assessment. Individual qualitative assessments of the 972 simulations were not performed. The relevance of this Thesis is in the implementation of the proposed algorithm to OCPs for which there are no known exact and analytical solutions. For these situations, the algorithm, with or without the aid of the ASAF structure, might provide a viable approach for obtaining acceptable candidate options to be recommended as a solution for the OCP2, saving considerable time and effort of the user by making efficient use of automation features. Therefore, the algorithm is presented as an important decision assistance tool for choosing an acceptable solution among many numeric possibilities found for a particular OCP. 2It is important to note that the algorithm does not guarantee optimal solution. In the case of the simulations present candidate solutions, the proposed algorithm does provide a recommended solution that could be acceptable or not, according to the decision of the user. Eventually, if an acceptable solution is not obtained in a number of contexts, the user may define new contexts to search for new solutions. 90 Chapter 5. DISCUSSION OF THE RESULTS The next chapter presents the conclusion of this doctoral Thesis, focusing on the importance of the obtained results. 91 6 CONCLUSION This Thesis presented the results of the research on the method to treat nonconvex optimal control problems by automated numerical simulations. After the study of the main references and theoretical concepts related to Optimization and Optimal Control, an algorithm to obtain numerical solutions recommended for nonconvex OCPs was proposed in Section 3.1. Furthermore, it was presented in Section 3.3, the ASAF structure, developed based on modern modeling, simulation and optimization tools, JModelica, (ANDERSSON, 2011; MAGNUSSON; ÅKESSON, 2012), IPOPT, (WÄCHTER; BIEGLER, 2006), and HSL routines, (Science and Technology Facilities Council, 2015), which was useful in the implementation of this algorithm. The verification of this Thesis hypothesis was performed through the study of four different nonconvex OCPs, in a total of 972 automated simulations. The resulting solutions were classified, verified and the application effectiveness was assessed, also, in an automated way. Thus, the algorithm was tested in a methodological manner using classical statistical methods that were able to estimate the median of its application effectiveness of 93%, with significance of 0.05 (confidence of 95%). Future works may explore the contexts that have already been defined in this research, performing additional simulations as an attempt to refute the hypothesis and, consequently, improve the proposed algorithm, or even develop new contexts for new OCPs of interest. 92 Chapter 6. CONCLUSION The research has impact when considering the range of OCPs for which it presents a viable approach to achieve solutions. The applications of optimal control theory in real physical problems, regardless if they are related to engineering, or biological systems, most frequently deal with the characteristics of non-compliance such as nonconvexity and nonlinearity, than with other problems that are already treated by the conventional theory. Therefore, with the results of this Thesis, there is an expansion of the possibilities for solving nonconvex OCP, which means an important and an original contribution to the subject. Other strong and interesting possibilities are related to the development of simulation, optimization and modeling tools, as well as computational techniques. Eventually, the application of the proposed method with the help of emerging technologies, when viable, may result in greater ease of the algorithm in obtaining acceptable recommended solutions. 93 REFERENCES ÅKESSON, J.; ÅRZÉN, K.-E.; GÄFVERT, M.; BERGDAHL, T.; TUMMESCHEIT, H. Modeling and Optimization with Optimica and JModelica.org—Languages and Tools for Solving Large-Scale Dynamic Optimization Problems. Computers and Chemical Engineering, v. 34, n. 11, p. 1737–1749, 2010. Available from Internet: . Cited on page 49. ANCONA, F.; BRESSAN, A.; CANNARSA, P.; CLARKE, F.; WOLENSKI, P. R. (Ed.). Geometric Control and Nonsmooth Analysis. In Honor of the 73rd Birthday of H. Hermes and of the 71st Birthday of R. T. Rockafellar. Proceedings of the conference, Rome, Italy, June 2006. [S.l.]: Hackensack, NJ: World Scientific, 2008. x + 361 p. ISBN 978-981-277-606-8/hbk. Cited 2 times on pages 25 and 32. ANDERSSON, J.; ÅKESSON, J.; CASELLA, F.; DIEHL, M. Integration of CasADi and JModelica.org. In: 8th International Modelica Conference 2011. Dresden, Germany: [s.n.], 2011. Available from Internet: . Cited 2 times on pages 49 and 91. ATHANS, M. Optimal control an introduction to the theory and its applications. New York: Dover Publications, 2007. ISBN 0486453286. Cited 2 times on pages 25 and 31. BANI-HANI, K.; GHABOUSSI, J. Nonlinear structural control using neural networks. Journal of Engineering Mechanics-asce, v. 124, n. 3, p. 319–327, mar. 1998. Cited on page 34. BAZARAA, M. S.; JARVIS, J. J.; SHERALI, H. D. Linear Programming and Network Flows. [S.l.]: Wiley, 2009. ISBN 0470462728. Cited on page 32. BERNSTEIN, D. Post-quantum Cryptography. Berlin: Springer, 2009. ISBN 9783642100192. Cited on page 88. BERTSEKAS, D. Dynamic programming and optimal control. Belmont, Mass: Athena Scientific, 2005. ISBN 1886529264. Cited on page 31. BERTSEKAS, D. Dynamic programming and optimal control. Belmont, Mass: Athena Scientific, 2005. ISBN 1886529264. Cited on page 31. BEWLEY, T. R.; TEMAN, R.; ZIANE, M. A general framework for robust control in fluid mechanics. Physica D, v. 138, n. 3-4, p. 360–392, abr. 2000. Cited on page 34. BONNARD, B.; CHYBA, M.; SUGNY, D. Time-minimal control of dissipative two-level quantum systems: The generic case. Ieee Transactions On Automatic Control, v. 54, n. 11, p. 2598–2610, nov. 2009. Cited on page 34. BRYSON, A. Applied optimal control. City: Hemisphere, 1975. ISBN 9780891162285. Cited 2 http://www.control.lth.se/documents/2010/ake+10cace.pdf http://www.control.lth.se/documents/2010/ake+10cace.pdf http://www.control.lth.se/documents/2011/and+11mod11.pdf 94 REFERENCES times on pages 25 and 31. CELLINA, A.; COLOMBO, G. On a Classical Problem of the Calculus of Variations without Convexity Assumptions. Annales de l’institut Henri Poincaré (C) Analyse non linéaire, Gauthier-Villars, v. 7, n. 2, p. 97–106, 1990. Cited 2 times on pages 25 and 32. CESARI, L. An Existence Theorem without Convexity Conditions. SIAM J. Control, Society for Industrial & Applied Mathematics, Philadelphia, v. 12, p. 319–331, 1974. ISSN 0036-1402. Cited 2 times on pages 25 and 32. CESARI, L. Optimization - Theory and Applications. Problems with Ordinary Differential Equations. 1983. Applications of Mathematics, 17. New York - Heidelberg - Berlin: Springer-Verlag. XIV, 542 p. Cited 3 times on pages 25, 32, and 42. CHESI, G.; GARULLI, A.; TESI, A.; VICINO, A. Solving quadratic distance problems: An lmi-based approach. Ieee Transactions On Automatic Control, v. 48, n. 2, p. 200–212, fev. 2003. Cited on page 34. CHOI, K. M.; CHO, S. W.; JUNG, H. J.; LEE, I. W. Semi-active fuzzy control for seismic response reduction using magnetorheological dampers. Earthquake Engineering & Structural Dynamics, v. 33, n. 6, p. 723–736, maio 2004. Cited on page 34. CIMEN, T.; BANKS, S. P. Global optimal feedback control for general nonlinear systems with nonquadratic performance criteria. Systems & Control Letters, v. 53, n. 5, p. 327–346, dez. 2004. Cited on page 34. CLARKE, F. A General Theorem on Necessary Conditions in Optimal Control. Discrete Contin. Dyn. Syst., American Institute of Mathematical Sciences (AIMS), Springfield, MO, v. 29, n. 2, p. 485–503, 2011. ISSN 1078-0947; 1553-5231/e. Cited 2 times on pages 25 and 32. CLARKE, F. H. Optimization and Nonsmooth Analysis. 1983. Canadian Mathematical Society Series of Monographs and Advanced Texts. A Wiley-Interscience Publication. New York: John Wiley & Sons, Inc. xiii, 308 pp. (1983). Cited 2 times on pages 25 and 32. DANIEL, W. Applied nonparametric statistics. Boston: PWS-KENT Pub, 1990. ISBN 0-534-91976-6. Cited 3 times on pages 29, 39, and 85. DING, Y.; JIA, Y.; WANG, S. S. Y. Identification of manning’s roughness coefficients in shallow water flows. Journal of Hydraulic Engineering-asce, v. 130, n. 6, p. 501–510, jun. 2004. Cited on page 34. ECO, U. Como se Faz Uma Tese (Em Português do Brasil). [S.l.]: Perspectiva, 2007. ISBN 8527300796. Cited 2 times on pages 31 and 87. EKELAND, I.; TEMAM, R. Convex Analysis and Variational Problems. Translated by Minerva Translations, Ltd., London. 1976. Studies in Mathematics and its Applications. Vol. 1. Amsterdam - Oxford: North-Holland Publishing Company; New York: American Elsevier Publishing Company, Inc. IX, 402 p. (1976). Cited 2 times on pages 25 and 32. ERKUS, B.; ABE, M.; FUJINO, Y. Investigation of semi-active control for seismic protection of elevated highway bridges. Engineering Structures, v. 24, n. 3, p. PII S0141–0296(01)00095–5, mar. 2002. Cited on page 34. REFERENCES 95 FRANKLIN, G. F.; ETC.; POWELL, D.; EMAMI-NAEINI, A. Feedback and Control of Dynamic Systems (Addison-Wesley series in electrical and computer engineering). 2nd. ed. [S.l.]: Addison Wesley Longman Publishing Co, 1991. ISBN 9780201508628. Cited on page 31. FRANKLIN, G. F.; POWELL, J. Digital Control of Dynamic Systems. [S.l.]: Addison Wesley, 1980. ISBN 9780201028911. Cited on page 31. GIBBONS, J. Nonparametric statistical inference. New York: Marcel Dekker, 2003. ISBN 0-8247-4052-1. Cited on page 39. GOTTWALD, G. A.; MELBOURNE, I. On the implementation of the 0–1 test for chaos. SIAM Journal on Applied Dynamical Systems, Society for Industrial & Applied Mathematics (SIAM), v. 8, n. 1, p. 129–145, jan 2009. Available from Internet: . Cited on page 55. HAGER, W. W. Runge-kutta methods in optimal control and the transformed adjoint system. Numerische Mathematik, v. 87, n. 2, p. 247–282, Dec 2000. ISSN 0945-3245. Available from Internet: . Cited on page 104. HOUSKA, B.; FERREAU, H. J.; DIEHL, M. An auto-generated real-time iteration algorithm for nonlinear mpc in the microsecond range. Automatica, v. 47, n. 10, p. 2279–2285, out. 2011. Cited on page 34. ILZHOEFER, A.; HOUSKA, B.; DIEHL, M. Nonlinear mpc of kites under varying wind conditions for a new class of large-scale wind power generators. International Journal of Robust and Nonlinear Control, v. 17, n. 17, p. 1590–1599, nov. 2007. Cited on page 34. JOO, Y. H.; SHIEH, L. S.; CHEN, G. R. Hybrid state-space fuzzy model-based controller with dual-rate sampling for digital control of chaotic systems. Ieee Transactions On Fuzzy Systems, v. 7, n. 4, p. 394–408, ago. 1999. Cited on page 34. KIRK, D. Optimal control theory; an introduction. Englewood Cliffs, N.J: Prentice-Hall, 1970. ISBN 0486434842. Cited 2 times on pages 25 and 31. KOCHE, J. C. Fundamentos de Metodologia Científica. Teoria da Ciência e Prática da Pesquisa. [S.l.]: Vozes, 2009. ISBN 8532618049. Cited on page 31. LENCI, S.; REGA, G. Optimal control of homoclinic bifurcation: Theoretical treatment and practical reduction of safe basin erosion in the helmholtz oscillator. Journal of Vibration and Control, v. 9, n. 3-4, p. 281–315, mar. 2003. Cited on page 34. LEVINE, D. M.; RAMSEY, P. P.; SMIDT, R. K. Applied Statistics for Engineers and Scientists: Using Microsoft Excel & Minitab. Pearson, 2000. ISBN 0134888014. Available from Internet: . Cited on page 38. LEWIS, F. L. Optimal Control. [S.l.]: Wiley-Interscience, 1986. ISBN 0471812404. Cited on page 31. LI, C.; FENG, G.; LIAO, X. Stabilization of nonlinear systems via periodically intermittent https://doi.org/10.1137/080718851 https://doi.org/10.1007/s002110000178 https://www.amazon.com/Applied-Statistics-Engineers-Scientists-Microsoft/dp/0134888014?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0134888014 https://www.amazon.com/Applied-Statistics-Engineers-Scientists-Microsoft/dp/0134888014?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0134888014 https://www.amazon.com/Applied-Statistics-Engineers-Scientists-Microsoft/dp/0134888014?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0134888014 96 REFERENCES control. Ieee Transactions On Circuits and Systems Ii-express Briefs, v. 54, n. 11, p. 1019–1023, nov. 2007. Cited on page 34. LIN, F.; ZHANG, W.; BRANDT, R. D. Robust hovering control of a pvtol aircraft. Ieee Transactions On Control Systems Technology, v. 7, n. 3, p. 343–351, maio 1999. Cited on page 34. LIU, H.; XIE, G.; WANG, L. Necessary and sufficient conditions for containment control of networked multi-agent systems. Automatica, v. 48, n. 7, p. 1415–1422, jul. 2012. Cited on page 34. LUENBERGER, D. G.; YE, Y. Linear and Nonlinear Programming (International Series in Operations Research & Management Science). Springer, 2008. ISBN 0387745025. Available from Internet: . Cited 2 times on pages 32 and 43. LUO, W. C.; CHU, Y. C.; LING, K. V. Inverse optimal adaptive control for attitude tracking of ’spacecraft. Ieee Transactions On Automatic Control, v. 50, n. 11, p. 1639–1654, nov. 2005. Cited on page 34. MACEK, W. M.; STRUMIK, M. Model for Hydromagnetic Convection in a Magnetized Fluid. Phys. Rev. E, American Physical Society, v. 82, n. 2, p. 27301, 2010. Available from Internet: . Cited on page 57. MACEK, W. M.; STRUMIK, M. Hyperchaotic Intermittent Convection in a Magnetized Viscous Fluid. Phys. Rev. Lett., American Physical Society, v. 112, n. 7, p. 74502, 2014. Available from Internet: . Cited on page 57. MAGNUSSON, F.; ÅKESSON, J. Collocation Methods for Optimization in a Modelica Environment. In: Proceedings of the 9th International Modelica Conference. [s.n.], 2012. Available from Internet: . Cited 11 times on pages 29, 32, 49, 91, 99, 100, 101, 102, 103, 104, and 105. MARCONI, M. d. A.; LAKATOS, E. M. Fundamentos de Metodologia Científica. [S.l.]: ATLAS - GRUPO GEN, 2010. ISBN 8522457581. Cited on page 31. MARICONDA, C. On a Parametric Problem of the Calculus of Variations without Convexity Assumptions. Journal of Mathematical Analysis and Applications, v. 170, n. 1, p. 291–297, 1992. ISSN 0022-247X. Available from Internet: . Cited 2 times on pages 25 and 32. MEZIAT, R.; PATIÑO, D.; PEDREGAL, P. An Alternative Approach for Non-linear Optimal Control Problems Based on the Method of Moments. Computational Optimization and Applications, v. 38, n. 1, p. 147–171, 2007. ISSN 1573-2894. Available from Internet: . Cited 3 times on pages 59, 62, and 64. OGATA, K. Discrete-Time Control Systems (2nd Edition). [S.l.]: Prentice Hall, 1995. ISBN 0133286428. Cited on page 31. OGATA, K. Engenharia de Controle Moderno. Pearson, 2010. ISBN 8587918230. Available https://www.amazon.com/Nonlinear-Programming-International-Operations-Management/dp/0387745025?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0387745025 https://www.amazon.com/Nonlinear-Programming-International-Operations-Management/dp/0387745025?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0387745025 https://www.amazon.com/Nonlinear-Programming-International-Operations-Management/dp/0387745025?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0387745025 https://www.amazon.com/Nonlinear-Programming-International-Operations-Management/dp/0387745025?SubscriptionId=0JYN1NVW651KCA56C102&tag=techkie-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=0387745025 https://link.aps.org/doi/10.1103/PhysRevE.82.027301 https://link.aps.org/doi/10.1103/PhysRevLett.112.074502 http://lup.lub.lu.se/record/2972280 http://www.sciencedirect.com/science/article/pii/0022247X9290020E http://www.sciencedirect.com/science/article/pii/0022247X9290020E http://dx.doi.org/10.1007/s10589-007-9032-1 http://dx.doi.org/10.1007/s10589-007-9032-1 REFERENCES 97 from Internet: . Cited on page 31. PARKER, T. S.; CHUA, L. Practical Numerical Algorithms for Chaotic Systems. [S.l.]: Springer, 1991. ISBN 0387966897. Cited 3 times on pages 35, 36, and 57. RAFIKOV, M.; BALTHAZAR, J. M. On control and synchronization in chaoti