Skip to content
Snippets Groups Projects
Commit d2e5fb85 authored by James Goppert's avatar James Goppert
Browse files

Simplex based trimming working through xml. 

See 737_cruise_steady_turn_simplex.xml script.
parent 5f3d1213
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
scicoslab/data/f16_lin.sce deleted 100644 → 0
+ 0
− 79
View file @ 5f3d1213
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);
src/FGFDMExec.cpp
+ 29
− 13
View file @ d2e5fb85
......@@ -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;
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......
src/Makefile.am
+ 1
− 19
View file @ d2e5fb85
......@@ -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 \
......
src/Trim.cpp deleted 100644 → 0
+ 0
− 318
View file @ 5f3d1213
/*
* 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
src/initialization/FGSimplexTrim.cpp 0 → 100644
+ 293
− 0
View file @ d2e5fb85
/*
* 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
src/initialization/FGSimplexTrim.h 0 → 100644
+ 81
− 0
View file @ d2e5fb85
/*
* 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
src/initialization/FGTrim.h
+ 2
− 1
View file @ d2e5fb85
......@@ -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>.
......
src/initialization/Makefile.am
+ 2
− 2
View file @ d2e5fb85
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment