diff --git a/scicoslab/data/f16_lin.sce b/scicoslab/data/f16_lin.sce deleted file mode 100644 index 99b2f626a7c1c5518c2e296b57fc0fa34b9453fb..0000000000000000000000000000000000000000 --- a/scicoslab/data/f16_lin.sce +++ /dev/null @@ -1,79 +0,0 @@ -f16.x0=.. -[ 5.0000021002e+02; - 3.9891191178e-02; - 3.9891191178e-02; - -2.3892269734e-05; - 3.9891191178e-02; - -2.3892269734e-05; - 7.7572299113e-03; - -6.9444185226e-18; - 1.8519409896e-07; - 7.3915343564e-09; - 1.0000000000e+03; - 1.5707963268e+00; - 0.0000000000e+00; - 0.0000000000e+00]; -f16.u0=.. -[ 2.6880598485e-01; - 6.7138558416e-03; - -2.6132360356e-02; - 1.4045552382e-01]; -f16.sys = syslin('c',.. -[ 1.3626752051e-01, 8.1722328147e+00, -3.2025072242e+01, -1.3351433275e+00, -3.2078892989e+01, -1.3351433275e+00, -1.0851775507e-01, -1.6335598151e-03, 1.7308220147e+00, 2.8257437912e-02, 4.5182250337e-05, -7.0213667056e-07, -1.0344734328e-11, 8.6483533293e-04; - -2.5722253396e-04, -9.6394940891e-01, 2.5278289277e-05, 9.9456442589e-01, -1.9383425850e-04, 9.9456442589e-01, -2.0174506277e-04, 6.8444220651e-08, 5.7073466616e-03, -1.9432890141e-04, 1.8968274778e-06, 1.1284964713e-06, 1.4456028966e-15, 1.2781709736e-08; - -4.7784519397e-08, -1.1293772630e-17, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 1.8534162488e-07, -1.4783071900e-08, 0.0000000000e+00, 0.0000000000e+00, 1.1417439440e-12, -1.1293772630e-17, 0.0000000000e+00, -1.0164395367e-16; - 6.2858088546e-05, -5.0565077129e+00, -4.5892127593e-03, 9.8681186356e-01, 3.5190006814e-02, 9.8681186356e-01, 2.3441127031e-02, -2.2686927199e-05, -2.6078069100e-01, -2.0592422564e-02, -1.1224312324e-05, 5.1895538648e-07, 0.0000000000e+00, -2.4957763756e-07; - -4.7784519397e-08, -1.1293772630e-17, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 1.8534162488e-07, -1.4783071900e-08, 0.0000000000e+00, 0.0000000000e+00, 1.1417439440e-12, -1.1293772630e-17, 0.0000000000e+00, -1.0164395367e-16; - 6.2858088546e-05, -5.0565077129e+00, -4.5892127593e-03, 9.8681186356e-01, 3.5190006814e-02, 9.8681186356e-01, 2.3441127031e-02, -2.2686927199e-05, -2.6078069100e-01, -2.0592422564e-02, -1.1224312324e-05, 5.1895538648e-07, 0.0000000000e+00, -2.4957763756e-07; - -3.2008894853e-06, -4.9396363888e-04, 4.9677102661e-04, -3.6327678409e-06, 4.9686121731e-04, -3.6327678409e-06, -3.2331999566e-01, 6.3846065510e-02, 4.0035218975e-02, -9.8827456506e-01, 1.5963872117e-06, 1.2629148306e-12, 0.0000000000e+00, 3.9202143757e-04; - 3.7097806918e-10, -4.4848987247e-13, 1.4806621203e-08, -2.6469779602e-19, 1.4806621203e-08, -2.6469779602e-19, 2.3911302360e-05, -1.9070788130e-06, 1.0000000000e+00, 3.9912364368e-02, -8.8636262811e-15, 9.7055858540e-19, 3.9704669403e-19, -6.6174449004e-19; - -6.1044155841e-05, -9.1442410075e-01, 1.3858413722e-06, -8.4986864823e-04, 6.4132592177e-07, -8.4986864823e-04, -2.8554997439e+01, 3.4354952575e-05, 3.3569450958e-01, 1.5481213535e+00, 2.1872259082e-04, -3.0195753305e-11, 0.0000000000e+00, -3.6667476997e-05; - 7.3173777182e-05, 2.1052609813e-02, 2.1757312387e-08, 2.8102302706e-04, -6.0211175397e-06, 2.8102302706e-04, 8.1164236010e+00, -1.7965656853e-05, -8.5077457351e-03, -8.0282553544e-01, -1.1372838184e-04, -5.9685994573e-09, 0.0000000000e+00, 1.9169024736e-05; - 1.5660216179e-11, -5.0000052090e+02, 4.9998537645e+02, 0.0000000000e+00, 4.9998537645e+02, 0.0000000000e+00, -5.3863163929e-12, 2.5957503891e-02, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, -1.9337829948e-11, 8.1926652962e-11, -1.5744928509e-10; - 1.4794832502e-11, -5.6318467818e-12, 1.8593212202e-07, -6.9469385111e-18, 1.8593212202e-07, -6.9469385111e-18, 9.5359737577e-07, -2.3947871511e-05, 0.0000000000e+00, 1.0007961815e+00, -3.5348598433e-16, 1.7090930090e-17, 3.7802153994e-18, -2.3892279770e-05; - 4.7784519371e-08, 2.5410988418e-17, 2.2587545260e-17, 0.0000000000e+00, 2.2587545260e-17, 0.0000000000e+00, -1.8534162480e-07, 7.3915375705e-09, 0.0000000000e+00, 0.0000000000e+00, -1.1417382971e-12, -1.8534162488e-07, 0.0000000000e+00, 1.1293772630e-17; - -3.7068293829e-10, 1.5440704768e-18, -6.2203982064e-18, 0.0000000000e+00, -6.2203982064e-18, 0.0000000000e+00, -2.3892279770e-05, 9.5283874322e-07, 0.0000000000e+00, 0.0000000000e+00, 8.8568103129e-15, -2.3892279770e-05, -1.5661286264e-18, 2.5146290622e-18],.. -[ 5.3830665042e+01, 7.6441978752e+01, 9.7117839711e+00, -5.4480264734e-04; - -4.2987608708e-03, 6.6914508611e-01, 3.9538168228e-02, -2.2179754683e-06; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - -1.7991826226e-01, 3.5451967364e-01, -7.1780342884e+00, 4.3702241171e-04; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - -1.7991826226e-01, 3.5451967364e-01, -7.1780342884e+00, 4.3702241171e-04; - -8.3517615734e-04, -9.1555664670e-03, -4.6330343352e-08, 1.2839806852e-03; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - -2.1134659127e-07, 1.6215240806e+01, 2.4327832709e-03, 1.8185195068e-01; - 6.6797898231e-08, 8.2722979134e-01, -1.0159303515e-04, -9.4311723130e-02; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00],.. -[ 9.9999999956e-01, 7.5791225148e-10, 3.7895612574e-10, 0.0000000000e+00, 3.7895612574e-10, 0.0000000000e+00, 1.8474111130e-09, 1.7053025658e-09, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 5.6843418861e-10, 0.0000000000e+00, 2.3684757859e-10; - -4.0476881106e-14, 1.0000000000e+00, -1.7925475918e-13, 0.0000000000e+00, -1.7925475918e-13, 0.0000000000e+00, 4.6259292693e-14, -3.9320398789e-13, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 9.8300996972e-14, 1.2721305490e-13, 4.0476881106e-14; - 0.0000000000e+00, -3.0289792254e-05, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 7.8031306637e-03, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 2.8912057933e-14, 3.0068540250e-13, -1.9081958236e-13; - -4.7784519397e-08, -1.1293772630e-17, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 1.8534162488e-07, -7.9056408410e-17, 0.0000000000e+00, 0.0000000000e+00, 1.1417439440e-12, -1.1293772630e-17, 0.0000000000e+00, -1.0164395367e-16; - 0.0000000000e+00, -3.0289792254e-05, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 7.8031306637e-03, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 2.8912057933e-14, 3.0068540250e-13, -1.9081958236e-13; - -4.7784519397e-08, -1.1293772630e-17, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 1.8534162488e-07, -7.9056408410e-17, 0.0000000000e+00, 0.0000000000e+00, 1.1417439440e-12, -1.1293772630e-17, 0.0000000000e+00, -1.0164395367e-16; - 1.4456028966e-15, 6.2883726004e-14, -9.0350181040e-14, 0.0000000000e+00, -9.0350181040e-14, 0.0000000000e+00, 1.0000000000e+00, -2.8044696195e-13, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, -3.4549909230e-13, 6.2160924556e-14, -2.0888961857e-13; - 0.0000000000e+00, 2.3148061657e-14, 5.0057736914e-13, 0.0000000000e+00, 5.0057736914e-13, 0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, -4.0798458820e-13, -1.0705978556e-13, 4.8900280430e-13; - 3.7038804218e-10, 2.2388810155e-13, -7.3915374614e-09, 0.0000000000e+00, -7.3915374614e-09, 0.0000000000e+00, 2.3873272313e-05, -9.5213839134e-07, 1.0000000000e+00, 0.0000000000e+00, -8.8496414142e-15, 1.7646519734e-19, 3.5293039469e-19, 9.5283874322e-07; - 1.4783062499e-11, -5.6094934734e-12, 1.8519417676e-07, 0.0000000000e+00, 1.8519417676e-07, 0.0000000000e+00, 9.5283874322e-07, 2.3855725072e-05, 0.0000000000e+00, 1.0000000000e+00, -3.5320474293e-16, -2.4167460230e-18, -1.3345180549e-18, -2.3873272313e-05; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.0000231365e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, -1.4802973662e-12, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.2952601954e-12, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 1.8503717077e-13, 1.4802973662e-12; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 1.0000000000e+00],.. -[ 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00; - 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00]); -f16.tfm = ss2tf(f16.sys); - diff --git a/src/FGFDMExec.cpp b/src/FGFDMExec.cpp index 59df63579481398ce27e270ab074fc4b2789f7aa..aa2870b28afbee4a6ea3317e11fb7c1d8e0054eb 100644 --- a/src/FGFDMExec.cpp +++ b/src/FGFDMExec.cpp @@ -62,6 +62,7 @@ INCLUDES //#include "initialization/FGTrimAnalysis.h" // Remove until later #include "input_output/FGPropertyManager.h" #include "input_output/FGScript.h" +#include "initialization/FGSimplexTrim.h" #include <iostream> #include <iterator> @@ -966,19 +967,34 @@ void FGFDMExec::DoTrim(int mode) void FGFDMExec::DoSimplexTrim(int mode) { - double saved_time; - - if (Constructing) return; - - if (mode < 0 || mode > JSBSim::tNone) { - cerr << endl << "Illegal trimming mode!" << endl << endl; - return; - } - saved_time = sim_time; - FGTrim trim(this, (JSBSim::TrimMode)mode); - if ( !trim.DoTrim() ) cerr << endl << "Trim Failed" << endl << endl; - trim.Report(); - sim_time = saved_time; + //double saved_time; + + //if (Constructing) return; + + //if (mode < 0 || mode > JSBSim::tNone) { + //cerr << endl << "Illegal trimming mode!" << endl << endl; + //return; + //} + //saved_time = sim_time; + //FGTrim trim(this, (JSBSim::TrimMode)mode); + //if ( !trim.DoTrim() ) cerr << endl << "Trim Failed" << endl << endl; + //trim.Report(); + //sim_time = saved_time; + double saved_time; + + if (Constructing) return; + + if (mode < 0 || mode > JSBSim::tNone) { + cerr << endl << "Illegal trimming mode!" << endl << endl; + return; + } + saved_time = sim_time; + + FGSimplexTrim trim(this, (JSBSim::TrimMode)mode); + //FGTrim trim(this, (JSBSim::TrimMode)mode); + //if ( !trim.DoTrim() ) cerr << endl << "Trim Failed" << endl << endl; + //trim.Report(); + sim_time = saved_time; } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/Makefile.am b/src/Makefile.am index e7a0c0ff47cef4334c538993f97ce6c4c67e459b..bfeba6861fbbb4c34dfc878ca1a84c10842147b8 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -9,7 +9,7 @@ LIBRARY_SOURCES = FGFDMExec.cpp FGJSBBase.cpp LIBRARY_INCLUDES = FGFDMExec.h FGJSBBase.h -noinst_PROGRAMS = JSBSim JSBSimTrim JSBSimFGComm +noinst_PROGRAMS = JSBSim JSBSimFGComm if BUILD_LIBRARIES @@ -34,9 +34,6 @@ libJSBSim_la_CXXFLAGS = $(AM_CXXFLAGS) JSBSim_SOURCES = JSBSim.cpp JSBSim_LDADD = libJSBSim.la -lm -JSBSimTrim_SOURCES = Trim.cpp -JSBSimTrim_LDADD = libJSBSim.la -lm - JSBSimFGComm_SOURCES = FGComm.cpp JSBSimFGComm_LDADD = libJSBSim.la -lm @@ -58,21 +55,6 @@ JSBSim_LDADD = \ simgear/magvar/libcoremag.a \ -lm -JSBSimTrim_SOURCES = Trim.cpp $(LIBRARY_SOURCES) -JSBSimTrim_LDADD = \ - initialization/libInit.a \ - models/atmosphere/libAtmosphere.a \ - models/libModels.a \ - models/flight_control/libFlightControl.a \ - models/propulsion/libPropulsion.a \ - input_output/libInputOutput.a \ - math/libMath.a \ - simgear/props/libProperties.a \ - simgear/xml/libExpat.a \ - simgear/magvar/libcoremag.a \ - -lm -Trim_CXXFLAGS=-Weffc++ - JSBSimFGComm_SOURCES = FGComm.cpp $(LIBRARY_SOURCES) JSBSimFGComm_LDADD = \ initialization/libInit.a \ diff --git a/src/Trim.cpp b/src/Trim.cpp deleted file mode 100644 index 6aa3fa17470702610856a695592d445c91829401..0000000000000000000000000000000000000000 --- a/src/Trim.cpp +++ /dev/null @@ -1,318 +0,0 @@ -/* - * Trim.cpp - * Copyright (C) James Goppert 2010 <james.goppert@gmail.com> - * - * Trim.cpp is free software: you can redistribute it and/or modify it - * under the terms of the GNU General Public License as published by the - * Free Software Foundation, either version 2 of the License, or - * (at your option) any later version. - * - * Trim.cpp is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - * See the GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License along - * with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include "initialization/FGTrimmer.h" -#include "math/FGStateSpace.h" -#include <iomanip> -#include <fstream> -#include "models/FGAircraft.h" -#include "models/propulsion/FGEngine.h" -#include "models/propulsion/FGTurbine.h" -#include "models/propulsion/FGTurboProp.h" -#include "math/FGNelderMead.h" -#include <stdexcept> -#include <fstream> - -template <class varType> -void prompt(const std::string & str, varType & var) -{ - std::cout << str + " [" << std::setw(10) << var << "]\t: "; - if (std::cin.peek() != '\n') - { - std::cin >> var; - std::cin.ignore(1000, '\n'); - } - else std::cin.get(); -} - -class Callback : public JSBSim::FGNelderMead::Callback -{ -private: - std::ofstream _outputFile; - JSBSim::FGTrimmer * _trimmer; -public: - Callback(std::string fileName, JSBSim::FGTrimmer * trimmer) : - _outputFile((fileName + std::string(".log")).c_str()), - _trimmer(trimmer) { - } - virtual ~Callback() { - _outputFile.close(); - } - void eval(const std::vector<double> &v) - { - _outputFile << _trimmer->eval(v) << std::endl;; - //std::cout << "v: "; - //for (int i=0;i<v.size();i++) std::cout << v[i] << " "; - //std::cout << std::endl; - } -}; - -int main (int argc, char const* argv[]) -{ - using namespace JSBSim; - - // variables - FGFDMExec fdm; - fdm.Setdt(1./120); - FGTrimmer::Constraints constraints; - - std::cout << "\n==============================================\n"; - std::cout << "\tJSBSim Trimming Utility\n"; - std::cout << "==============================================\n" << std::endl; - - // defaults - constraints.velocity = 500; - std::string aircraft="f16"; - double rtol = std::numeric_limits<float>::epsilon(); - double abstol = 1e-2;//std::numeric_limits<double>::epsilon(); - double speed = 1.1; // > 1 - double random = 0; // random scale factor added to all simplex calcs - int iterMax = 2000; - bool showConvergeStatus = false; - bool pause = false; - bool showSimplex = false; - bool variablePropPitch = false; - int debugLevel = 0; - std::string fileName = aircraft; - - // input - std::cout << "input ( press enter to accept [default] )\n" << std::endl; - - // load model - std::string aircraftName = ""; - prompt("\tdebug level\t\t",debugLevel); - fdm.SetDebugLevel(debugLevel); - std::cout << "model selection" << std::endl; - while (1) - { - prompt("\taircraft\t\t",aircraft); - prompt("\toutput file name\t",fileName); - fdm.LoadModel("../aircraft","../engine","../systems",aircraft); - aircraftName = fdm.GetAircraft()->GetAircraftName(); - if (aircraftName == "") - { - std::cout << "\tfailed to load aircraft" << std::endl; - } - else - { - std::cout << "\tsuccessfully loaded: " << aircraftName << std::endl; - break; - } - } - - // Turn on propulsion system - fdm.GetPropulsion()->InitRunning(-1); - - // get propulsion pointer to determine type/ etc. - FGEngine * engine0 = fdm.GetPropulsion()->GetEngine(0); - FGThruster * thruster0 = engine0->GetThruster(); - - // flight conditions - std::cout << "\nflight conditions: " << std::endl; - prompt("\taltitude, ft\t\t",constraints.altitude); - prompt("\tvelocity, ft/s\t\t",constraints.velocity); - prompt("\tgamma, deg\t\t",constraints.gamma); - if (thruster0->GetType()==FGThruster::ttPropeller) - prompt("\tvariable prop pitch?\t\t",variablePropPitch); - constraints.gamma *= M_PI/180; - - // mode menu - while (1) - { - int mode = 0; - prompt("\tmode < non-turning(0), rolling(1), pitching(2), yawing(3) >",mode); - constraints.rollRate = 0; - constraints.pitchRate = 0; - constraints.yawRate = 0; - if (mode == 0) break; - else if (mode == 1) - { - prompt("\troll rate, rad/s",constraints.rollRate); - prompt("\tstability axis roll",constraints.stabAxisRoll); - break; - } - else if (mode == 2) - { - prompt("\tpitch rate, rad/s",constraints.pitchRate); - break; - } - else if (mode == 3) - { - prompt("\tyaw rate, rad/s",constraints.yawRate); - break; - } - else std::cout << "\tunknown mode: " << mode << std::endl; - } - - // solver properties - std::cout << "\nsolver properties: " << std::endl; - std::cout << std::scientific; - prompt("\tshow converge status?\t",showConvergeStatus); - prompt("\tshow simplex?\t\t",showSimplex); - prompt("\tpause?\t\t\t",pause); - prompt("\trelative tolerance\t",rtol); - prompt("\tabsolute tolerance\t",abstol); - prompt("\tmax iterations\t\t",iterMax); - prompt("\tconvergence speed\t",speed); - prompt("\trandomization ratio\t",random); - std::cout << std::fixed; - - // initial solver state - int n = 6; - std::vector<double> initialGuess(n), lowerBound(n), upperBound(n), initialStepSize(n); - - lowerBound[0] = 0; //throttle - lowerBound[1] = -1; // elevator - lowerBound[2] = -90*M_PI/180; // alpha - lowerBound[3] = -1; // aileron - lowerBound[4] = -1; // rudder - lowerBound[5] = -90*M_PI/180; // beta - - upperBound[0] = 1; //throttle - upperBound[1] = 1; // elevator - upperBound[2] = 90*M_PI/180; // alpha - upperBound[3] = 1; // aileron - upperBound[4] = 1; // rudder - upperBound[5] = 90*M_PI/180; // beta - - initialStepSize[0] = 0.2; //throttle - initialStepSize[1] = 0.1; // elevator - initialStepSize[2] = 0.1; // alpha - initialStepSize[3] = 0.1; // aileron - initialStepSize[4] = 0.1; // rudder - initialStepSize[5] = 0.1; // beta - - initialGuess[0] = 0.5; // throttle - initialGuess[1] = 0; // elevator - initialGuess[2] = 0; // alpha - initialGuess[3] = 0; // aileron - initialGuess[4] = 0; // rudder - initialGuess[5] = 0; // beta - - // solve - FGTrimmer trimmer(fdm, constraints); - Callback callback(fileName,&trimmer); - FGNelderMead * solver; - try - { - solver = new FGNelderMead(trimmer,initialGuess, - lowerBound, upperBound, initialStepSize,iterMax,rtol, - abstol,speed,random,showConvergeStatus,showSimplex,pause,&callback); - while(solver->status()==1) solver->update(); - } - catch (const std::runtime_error & e) - { - std::cout << e.what() << std::endl; - return 1; - } - - // output - try - { - trimmer.printSolution(solver->getSolution()); // this also loads the solution into the fdm - } - catch(std::runtime_error & e) - { - std::cout << "caught std::runtime error" << std::endl; - std::cout << "exception: " << e.what() << std::endl; - return 1; - } - - //std::cout << "\nsimulating flight to determine trim stability" << std::endl; - - //std::cout << "\nt = 5 seconds" << std::endl; - //for (int i=0;i<5*120;i++) fdm.Run(); - //trimmer.printState(); - - //std::cout << "\nt = 10 seconds" << std::endl; - //for (int i=0;i<5*120;i++) fdm.Run(); - //trimmer.printState(); - - std::cout << "\nlinearization: " << std::endl; - FGStateSpace ss(fdm); - - ss.x.add(new FGStateSpace::Vt); - ss.x.add(new FGStateSpace::Alpha); - ss.x.add(new FGStateSpace::Theta); - ss.x.add(new FGStateSpace::Q); - ss.x.add(new FGStateSpace::Theta); - ss.x.add(new FGStateSpace::Q); - - if (thruster0->GetType()==FGThruster::ttPropeller) - { - ss.x.add(new FGStateSpace::Rpm0); - if (variablePropPitch) ss.x.add(new FGStateSpace::PropPitch); - int numEngines = fdm.GetPropulsion()->GetNumEngines(); - if (numEngines>1) ss.x.add(new FGStateSpace::Rpm1); - if (numEngines>2) ss.x.add(new FGStateSpace::Rpm2); - if (numEngines>3) ss.x.add(new FGStateSpace::Rpm3); - } - ss.x.add(new FGStateSpace::Beta); - ss.x.add(new FGStateSpace::Phi); - ss.x.add(new FGStateSpace::P); - ss.x.add(new FGStateSpace::R); - ss.x.add(new FGStateSpace::Alt); - ss.x.add(new FGStateSpace::Psi); - ss.x.add(new FGStateSpace::Longitude); - ss.x.add(new FGStateSpace::Latitude); - - ss.u.add(new FGStateSpace::ThrottleCmd); - ss.u.add(new FGStateSpace::DaCmd); - ss.u.add(new FGStateSpace::DeCmd); - ss.u.add(new FGStateSpace::DrCmd); - - // state feedback - ss.y = ss.x; - - std::vector< std::vector<double> > A,B,C,D; - std::vector<double> x0 = ss.x.get(), u0 = ss.u.get(); - std::vector<double> y0 = x0; // state feedback - std::cout << ss << std::endl; - - ss.linearize(x0,u0,y0,A,B,C,D); - - int width=10; - std::cout.precision(3); - std::cout - << std::fixed - << std::right - << "\nA=\n" << std::setw(width) << A - << "\nB=\n" << std::setw(width) << B - << "\nC=\n" << std::setw(width) << C - << "\nD=\n" << std::setw(width) << D - << std::endl; - - // write scicoslab file - std::ofstream scicos(std::string(aircraft+"_lin.sce").c_str()); - scicos.precision(10); - width=20; - scicos - << std::scientific - << aircraft << ".x0=..\n" << std::setw(width) << x0 << ";\n" - << aircraft << ".u0=..\n" << std::setw(width) << u0 << ";\n" - << aircraft << ".sys = syslin('c',..\n" - << std::setw(width) << A << ",..\n" - << std::setw(width) << B << ",..\n" - << std::setw(width) << C << ",..\n" - << std::setw(width) << D << ");\n" - << aircraft << ".tfm = ss2tf(" << aircraft << ".sys);\n" - << std::endl; -} - - -// vim:ts=4:sw=4 diff --git a/src/initialization/FGSimplexTrim.cpp b/src/initialization/FGSimplexTrim.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d2fe7c150132ede0a96b151fe390167e5ce8b3b1 --- /dev/null +++ b/src/initialization/FGSimplexTrim.cpp @@ -0,0 +1,293 @@ +/* + * FGSimplexTrim.cpp + * Copyright (C) James Goppert 2010 <james.goppert@gmail.com> + * + * FGSimplexTrim.cpp is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * FGSimplexTrim.cpp is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "FGSimplexTrim.h" + +namespace JSBSim { + +FGSimplexTrim::FGSimplexTrim(FGFDMExec * fdmPtr, TrimMode mode) +{ + using namespace JSBSim; + + // variables + FGFDMExec & fdm = *fdmPtr; + fdm.Setdt(1./120); + FGTrimmer::Constraints constraints; + + std::cout << "\n==============================================\n"; + std::cout << "\tJSBSim Trimming Utility\n"; + std::cout << "==============================================\n" << std::endl; + + // defaults + constraints.velocity = fdm.GetAuxiliary()->GetVcalibratedFPS(); + constraints.altitude = fdm.GetPropagate()->GetAltitudeASL(); + std::string aircraft = fdm.GetAircraft()->GetAircraftName(); + double rtol = std::numeric_limits<float>::epsilon(); + double abstol = 1e-2;//std::numeric_limits<double>::epsilon(); + double speed = 1.1; // > 1 + double random = 0; // random scale factor added to all simplex calcs + int iterMax = 2000; + bool showConvergeStatus = false; + bool pause = false; + bool showSimplex = false; + bool variablePropPitch = false; + int debugLevel = 0; + std::string fileName = aircraft; + + // input + std::cout << "input ( press enter to accept [default] )\n" << std::endl; + + // load model + std::string aircraftName = fdm.GetAircraft()->GetAircraftName(); + //prompt("\tdebug level\t\t",debugLevel); + //fdm.SetDebugLevel(debugLevel); + //std::cout << "model selection" << std::endl; + //while (1) + //{ + //prompt("\taircraft\t\t",aircraft); + //prompt("\toutput file name\t",fileName); + //fdm.LoadModel("../aircraft","../engine","../systems",aircraft); + //aircraftName = fdm.GetAircraft()->GetAircraftName(); + //if (aircraftName == "") + //{ + //std::cout << "\tfailed to load aircraft" << std::endl; + //} + //else + //{ + //std::cout << "\tsuccessfully loaded: " << aircraftName << std::endl; + //break; + //} + //} + + // Turn on propulsion system + fdm.GetPropulsion()->InitRunning(-1); + + // get propulsion pointer to determine type/ etc. + FGEngine * engine0 = fdm.GetPropulsion()->GetEngine(0); + FGThruster * thruster0 = engine0->GetThruster(); + + // flight conditions + //std::cout << "\nflight conditions: " << std::endl; + //prompt("\taltitude, ft\t\t",constraints.altitude); + //prompt("\tvelocity, ft/s\t\t",constraints.velocity); + //prompt("\tgamma, deg\t\t",constraints.gamma); + + // TODO check that this works properly + constraints.gamma = fdm.GetPropagate()->GetEuler(2); + + //if (thruster0->GetType()==FGThruster::ttPropeller) + //prompt("\tvariable prop pitch?\t\t",variablePropPitch); + // FIXME, enable + constraints.gamma *= M_PI/180; + + // mode menu + while (1) + { + //prompt("\tmode",mode); + constraints.rollRate = 0; + constraints.pitchRate = 0; + constraints.yawRate = 0; + if (mode == tLongitudinal) break; + else if (mode == tRoll) + { + //prompt("\troll rate, rad/s",constraints.rollRate); + //prompt("\tstability axis roll",constraints.stabAxisRoll); + // TODO check that this works properly + constraints.pitchRate = fdm.GetPropagate()->GetPQR(1); + constraints.stabAxisRoll = true; // FIXME, make this an option + break; + } + else if (mode == tPullup) + { + //prompt("\tpitch rate, rad/s",constraints.pitchRate); + // TODO check that this works properly + constraints.pitchRate = fdm.GetPropagate()->GetPQR(2); + break; + } + else if (mode == tTurn) + { + //prompt("\tyaw rate, rad/s",constraints.yawRate); + // TODO check that this works properly + constraints.yawRate = fdm.GetPropagate()->GetPQR(3); + break; + } + else { + std::cerr << "\tunknown mode: " << mode << std::endl; + exit(1); + } + } + + // solver properties + // TODO make these options + //std::cout << "\nsolver properties: " << std::endl; + //std::cout << std::scientific; + //prompt("\tshow converge status?\t",showConvergeStatus); + //prompt("\tshow simplex?\t\t",showSimplex); + //prompt("\tpause?\t\t\t",pause); + //prompt("\trelative tolerance\t",rtol); + //prompt("\tabsolute tolerance\t",abstol); + //prompt("\tmax iterations\t\t",iterMax); + //prompt("\tconvergence speed\t",speed); + //prompt("\trandomization ratio\t",random); + //std::cout << std::fixed; + + // initial solver state + int n = 6; + std::vector<double> initialGuess(n), lowerBound(n), upperBound(n), initialStepSize(n); + + lowerBound[0] = 0; //throttle + lowerBound[1] = -1; // elevator + lowerBound[2] = -90*M_PI/180; // alpha + lowerBound[3] = -1; // aileron + lowerBound[4] = -1; // rudder + lowerBound[5] = -90*M_PI/180; // beta + + upperBound[0] = 1; //throttle + upperBound[1] = 1; // elevator + upperBound[2] = 90*M_PI/180; // alpha + upperBound[3] = 1; // aileron + upperBound[4] = 1; // rudder + upperBound[5] = 90*M_PI/180; // beta + + initialStepSize[0] = 0.2; //throttle + initialStepSize[1] = 0.1; // elevator + initialStepSize[2] = 0.1; // alpha + initialStepSize[3] = 0.1; // aileron + initialStepSize[4] = 0.1; // rudder + initialStepSize[5] = 0.1; // beta + + initialGuess[0] = 0.5; // throttle + initialGuess[1] = 0; // elevator + initialGuess[2] = 0; // alpha + initialGuess[3] = 0; // aileron + initialGuess[4] = 0; // rudder + initialGuess[5] = 0; // beta + + // solve + FGTrimmer trimmer(fdm, constraints); + Callback callback(fileName,&trimmer); + FGNelderMead * solver; + try + { + solver = new FGNelderMead(trimmer,initialGuess, + lowerBound, upperBound, initialStepSize,iterMax,rtol, + abstol,speed,random,showConvergeStatus,showSimplex,pause,&callback); + while(solver->status()==1) solver->update(); + } + catch (const std::runtime_error & e) + { + std::cout << e.what() << std::endl; + exit(1); + } + + // output + try + { + trimmer.printSolution(solver->getSolution()); // this also loads the solution into the fdm + } + catch(std::runtime_error & e) + { + std::cout << "caught std::runtime error" << std::endl; + std::cout << "exception: " << e.what() << std::endl; + exit(1); + } + + //std::cout << "\nsimulating flight to determine trim stability" << std::endl; + + //std::cout << "\nt = 5 seconds" << std::endl; + //for (int i=0;i<5*120;i++) fdm.Run(); + //trimmer.printState(); + + //std::cout << "\nt = 10 seconds" << std::endl; + //for (int i=0;i<5*120;i++) fdm.Run(); + //trimmer.printState(); + + std::cout << "\nlinearization: " << std::endl; + FGStateSpace ss(fdm); + + ss.x.add(new FGStateSpace::Vt); + ss.x.add(new FGStateSpace::Alpha); + ss.x.add(new FGStateSpace::Theta); + ss.x.add(new FGStateSpace::Q); + ss.x.add(new FGStateSpace::Theta); + ss.x.add(new FGStateSpace::Q); + + if (thruster0->GetType()==FGThruster::ttPropeller) + { + ss.x.add(new FGStateSpace::Rpm0); + if (variablePropPitch) ss.x.add(new FGStateSpace::PropPitch); + int numEngines = fdm.GetPropulsion()->GetNumEngines(); + if (numEngines>1) ss.x.add(new FGStateSpace::Rpm1); + if (numEngines>2) ss.x.add(new FGStateSpace::Rpm2); + if (numEngines>3) ss.x.add(new FGStateSpace::Rpm3); + } + ss.x.add(new FGStateSpace::Beta); + ss.x.add(new FGStateSpace::Phi); + ss.x.add(new FGStateSpace::P); + ss.x.add(new FGStateSpace::R); + ss.x.add(new FGStateSpace::Alt); + ss.x.add(new FGStateSpace::Psi); + ss.x.add(new FGStateSpace::Longitude); + ss.x.add(new FGStateSpace::Latitude); + + ss.u.add(new FGStateSpace::ThrottleCmd); + ss.u.add(new FGStateSpace::DaCmd); + ss.u.add(new FGStateSpace::DeCmd); + ss.u.add(new FGStateSpace::DrCmd); + + // state feedback + ss.y = ss.x; + + std::vector< std::vector<double> > A,B,C,D; + std::vector<double> x0 = ss.x.get(), u0 = ss.u.get(); + std::vector<double> y0 = x0; // state feedback + std::cout << ss << std::endl; + + ss.linearize(x0,u0,y0,A,B,C,D); + + int width=10; + std::cout.precision(3); + std::cout + << std::fixed + << std::right + << "\nA=\n" << std::setw(width) << A + << "\nB=\n" << std::setw(width) << B + << "\nC=\n" << std::setw(width) << C + << "\nD=\n" << std::setw(width) << D + << std::endl; + + // write scicoslab file + std::ofstream scicos(std::string(aircraft+"_lin.sce").c_str()); + scicos.precision(10); + width=20; + scicos + << std::scientific + << aircraft << ".x0=..\n" << std::setw(width) << x0 << ";\n" + << aircraft << ".u0=..\n" << std::setw(width) << u0 << ";\n" + << aircraft << ".sys = syslin('c',..\n" + << std::setw(width) << A << ",..\n" + << std::setw(width) << B << ",..\n" + << std::setw(width) << C << ",..\n" + << std::setw(width) << D << ");\n" + << aircraft << ".tfm = ss2tf(" << aircraft << ".sys);\n" + << std::endl; +} + +} // JSBSim + +// vim:ts=4:sw=4 diff --git a/src/initialization/FGSimplexTrim.h b/src/initialization/FGSimplexTrim.h new file mode 100644 index 0000000000000000000000000000000000000000..31c14daf0a5ea3e67ea65beffeb8314d497d654c --- /dev/null +++ b/src/initialization/FGSimplexTrim.h @@ -0,0 +1,81 @@ +/* + * FGSimplexTrim.h + * Copyright (C) James Goppert 2010 <james.goppert@gmail.com> + * + * FGSimplexTrim.h is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * FGSimplexTrim.h is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef FGSimplexTrim_H_ +#define FGSimplexTrim_H_ + +#include "initialization/FGTrimmer.h" +#include "math/FGStateSpace.h" +#include <iomanip> +#include <fstream> +#include "models/FGAircraft.h" +#include "models/propulsion/FGEngine.h" +#include "models/propulsion/FGTurbine.h" +#include "models/propulsion/FGTurboProp.h" +#include "math/FGNelderMead.h" +#include <stdexcept> +#include <fstream> +#include <cstdlib> + +namespace JSBSim { + +class FGSimplexTrim +{ +public: + FGSimplexTrim(FGFDMExec * fdmPtr, TrimMode mode); +private: + template <class varType> + void prompt(const std::string & str, varType & var) + { + std::cout << str + " [" << std::setw(10) << var << "]\t: "; + if (std::cin.peek() != '\n') + { + std::cin >> var; + std::cin.ignore(1000, '\n'); + } + else std::cin.get(); + } + + class Callback : public JSBSim::FGNelderMead::Callback + { + private: + std::ofstream _outputFile; + JSBSim::FGTrimmer * _trimmer; + public: + Callback(std::string fileName, JSBSim::FGTrimmer * trimmer) : + _outputFile((fileName + std::string(".log")).c_str()), + _trimmer(trimmer) { + } + virtual ~Callback() { + _outputFile.close(); + } + void eval(const std::vector<double> &v) + { + _outputFile << _trimmer->eval(v) << std::endl;; + //std::cout << "v: "; + //for (int i=0;i<v.size();i++) std::cout << v[i] << " "; + //std::cout << std::endl; + } + }; +}; + +} // JSBSim + +#endif //FGSimplexTrim_H_ + +// vim:ts=4:sw=4 diff --git a/src/initialization/FGTrim.h b/src/initialization/FGTrim.h index 0d6e3092cdd260745df07ada44d3edfa6082bf19..bfab7ccf2787116f796521a104d979bde5b71eda 100644 --- a/src/initialization/FGTrim.h +++ b/src/initialization/FGTrim.h @@ -69,7 +69,7 @@ FORWARD DECLARATIONS namespace JSBSim { typedef enum { tLongitudinal=0, tFull, tGround, tPullup, - tCustom, tTurn, tNone } TrimMode; + tCustom, tTurn, tRoll, tNone } TrimMode; /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CLASS DOCUMENTATION @@ -94,6 +94,7 @@ CLASS DOCUMENTATION - tPullup: tLongitudinal but adjust alpha to achieve load factor input with SetTargetNlf() - tGround: wdot with altitude, qdot with theta, and pdot with phi + - tRoll: trim to a constant roll rate The remaining modes include <b>tCustom</b>, which is completely user defined and <b>tNone</b>. diff --git a/src/initialization/Makefile.am b/src/initialization/Makefile.am index 1e908e60c2b986aa1c866d84b336be9b52226d95..6fdf6fa0aba60b7393dffb335890f494b78bbde2 100644 --- a/src/initialization/Makefile.am +++ b/src/initialization/Makefile.am @@ -2,10 +2,10 @@ includedir = @includedir@/JSBSim/initialization ###AM_CPPFLAGS = -DOLD_LIBC -DAGO_DIRECTSEARCH -Wno-non-template-friend -LIBRARY_SOURCES = FGInitialCondition.cpp FGTrim.cpp FGTrimAxis.cpp FGTrimmer.cpp +LIBRARY_SOURCES = FGInitialCondition.cpp FGTrim.cpp FGTrimAxis.cpp FGTrimmer.cpp FGSimplexTrim.cpp ### FGTrimAnalysis.cpp FGTrimAnalysisControl.cpp -LIBRARY_INCLUDES = FGInitialCondition.h FGTrim.h FGTrimAxis.h FGTrimmer.h +LIBRARY_INCLUDES = FGInitialCondition.h FGTrim.h FGTrimAxis.h FGTrimmer.h FGSimplexTrim.h ### FGTrimAnalysis.h FGTrimAnalysisControl.h if BUILD_LIBRARIES