PS_Method_CCS_test/Pseudo_Spectral_Method/Orion_example/OrionShipTest.cpp

#include <iostream>
#include <armadillo>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include "GPM.h"

// 用户新增头文件
#include "OrionShip.h"

using namespace std;
using namespace arma;

int main_Orion() {
	// 计时器开始
	double dur;
	clock_t start, end;
	start = clock();

	////////////////////////////////////参数初始化//////////////////////////////////////////////
	input0 input;
	output0* output;
	auxdata1 auxdata;
	// 伪谱法条件初始化
	Intput_Initial(&input);

	// 人为编辑定义部分
	auxdata.Sref = 23.8;
	auxdata.M0 = 9500;
	auxdata.C1_v = 2e-8;		// 与热流密度相关参数
	auxdata.Rd_v = 0.006;		// 与热流密度相关参数
	auxdata.n_v = 0.5;			// 与热流密度相关参数
	auxdata.m_v = 3.15;
	auxdata.g = 9.80665;
	auxdata.Re = 6371110;
	input.setup.auxdata = auxdata;

	// 归一化单位向量
	double pi = datum::pi;
	double g = auxdata.g;
	double Re = auxdata.Re;
	double tc = sqrt(Re / g);
	double rc = Re;
	double vc = sqrt(Re * g);
	vec state_c; state_c << vc << 1.0 << 1.0 << rc << 1.0 << 1.0 << endr;

	// 范围赋值
	int length_state = 6;
	int length_control = 1;

	double initialtime_lower = 0;
	double initialtime_upper = 0;

	double finaltime_lower = 0;
	double finaltime_upper = 4500.0 / tc;

	vec initialstate_lower; initialstate_lower << 10600.0 << -5.9 / 180.0 * pi << 90.0 / 180.0 * pi  << 120.0 * 1e3 + Re << 0.0 / 180.0 * pi   << 0.0 / 180.0 * pi << endr;
	vec initialstate_upper; initialstate_upper << 10600.0 << -5.9 / 180.0 * pi << 90.0 / 180.0 * pi  << 120.0 * 1e3 + Re << 0.0 / 180.0 * pi   << 0.0 / 180.0 * pi << endr;

	vec state_lower; state_lower              <<   100.0 << -80.0 / 180.0 * pi << -90.0 / 180.0 * pi <<  10.0 * 1e3 + Re <<  -1.0 / 180.0 * pi << -90.0 / 180.0 * pi << endr;
	vec state_upper; state_upper              << 11700.0 << +20.0 / 180.0 * pi << 270.0 / 180.0 * pi << 200.0 * 1e3 + Re << 180.0 / 180.0 * pi << +90.0 / 180.0 * pi << endr;

	vec finalstate_lower; finalstate_lower   << 100.0    << -80.0 / 180.0 * pi << -90.0 / 180.0 * pi << 10.0 * 1e3 + Re  << 90.0 / 180.0 * pi  << -90.0 / 180.0 * pi << endr;
	vec finalstate_upper; finalstate_upper   << 300.0    << + 0.0 / 180.0 * pi << 270.0 / 180.0 * pi << 10.0 * 1e3 + Re  << 90.0 / 180.0 * pi  << +90.0 / 180.0 * pi << endr;

	vec control_lower; control_lower << -180.0 / 180.0 * pi << endr;
	vec control_upper; control_upper << +180.0 / 180.0 * pi << endr;

	vec path_lower; path_lower << -6.5 <<  endr;
	vec path_upper; path_upper << 6.5 << endr;

	rowvec integral_lower; integral_lower.reset();
	rowvec integral_upper; integral_upper.reset();

	input.setup.state_length = length_state;
	input.setup.control_length = length_control;
	input.setup.bounds.phase.initialtime.lower = initialtime_lower;
	input.setup.bounds.phase.initialtime.upper = initialtime_upper;
	input.setup.bounds.phase.finaltime.lower = finaltime_lower;
	input.setup.bounds.phase.finaltime.upper = finaltime_upper;
	input.setup.bounds.phase.initialstate.lower = initialstate_lower / state_c;
	input.setup.bounds.phase.initialstate.upper = initialstate_upper / state_c;
	input.setup.bounds.phase.state.lower = state_lower / state_c;
	input.setup.bounds.phase.state.upper = state_upper / state_c;        	// 默认状态上限为正无穷
	input.setup.bounds.phase.finalstate.lower = finalstate_lower / state_c;
	input.setup.bounds.phase.finalstate.upper = finalstate_upper / state_c;
	input.setup.bounds.phase.control.lower = control_lower;
	input.setup.bounds.phase.control.upper = control_upper;
	input.setup.bounds.phase.path.lower = path_lower;
	input.setup.bounds.phase.path.upper = path_upper;
	input.setup.bounds.phase.integral.lower = integral_lower;
	input.setup.bounds.phase.integral.upper = integral_upper;
	input.setup.bounds.phase.event.lower.reset();
	input.setup.bounds.phase.event.upper.reset();

	// 猜测值
	guess0 guess;
	int length_path = 1;
	vec guess_time; guess_time << 0.0 << 2500.0 / tc << endr;
	mat guess_state(2, length_state, fill::zeros);
	mat guess_control(2, length_control, fill::zeros);
	mat guess_path(2, length_path, fill::zeros);
	guess_state.row(0) = input.setup.bounds.phase.initialstate.lower.t();
	guess_state.row(1) = input.setup.bounds.phase.finalstate.lower.t();
	guess_control.row(0) = input.setup.bounds.phase.control.lower.t();
	guess_control.row(1) = input.setup.bounds.phase.control.upper.t();
	guess_path.row(0) = input.setup.bounds.phase.path.lower.t();
	guess_path.row(1) = input.setup.bounds.phase.path.upper.t();

	guess.phase.time = guess_time;
	guess.phase.state = guess_state;
	guess.phase.control = guess_control;
	guess.phase.path = guess_path;
	input.setup.guess = guess;

	// 配点设置
	int inter_num = 5;
	input.setup.name = "OrionShip";							// 问题
	mesh1 mesh;
	mesh.method = "p-h";								// 网格细化方法的名字，字符串表示
	mesh.tolerance = 1e-5;								// 网格细化精度的允许值
	mesh.maxiterration = 6;								// 网格细化最大迭代次数
	mesh.colpointsmin = 2;								// 网格细化产生新interval下的最小配点数目
	mesh.colpointsmax = 20;								// 网格细化产生新interval下的最大配点数目
	mesh.splitmult = 1.2;								// 比例因子，确定创建新网格的间隔数目
	mesh.curveratio = 2;								// 网格间隔中解的最大平均曲率比的阈值，以确定是否增加新间隔
	mesh.phase.inter_num = inter_num;					// 段p内的间隔数目, 整数
	mesh.phase.colpoints = 2 * ones(inter_num, 1);		// 段p内的配置点数目，整数数组，长度与间隔数目相等
	mesh.phase.i_proportion = 2.0 / inter_num * ones(inter_num, 1); // 初始网格划分的比例（当前为均等划分，和为2）
	input.setup.mesh = mesh;							// 组合setup.mesh结构体
	// input.setup.Scal = 0;
	input.setup.warmstart = 0;
	input.setup.interp.method = "linear";

	// 函数设置 
	functions1 functions;
	functions.continous1 = Reentry_Orion_Dynamic;
	functions.endpoint = OrionObjective;

	// 求解器设置
	NLP1 NLP;
	NLP.ipoptoptions.print_level = 1;
	NLP.ipoptoptions.hessian_approximation = "limited-memory";
	NLP.ipoptoptions.mu_strategy = "adaptive";
	NLP.ipoptoptions.linear_solver = "ma27"; //mumps、pardiso
	NLP.ipoptoptions.tolerance = 1e-1;
	NLP.ipoptoptions.maxiterration = 3000;
	input.setup.NLP = NLP;
	input.setup.derivatives.supplier = "CD";

	//input.setup.feature.FSE = 1;	///////

	// 伪谱法过程实现
	int iteration_final;
	output = PM_Process(&input, functions, &iteration_final);

	// 计时结束
	end = clock();
	dur = (double)(end - start);
#ifndef LUOYC20220721
	printf("Use Time: %f \n", (dur / CLOCKS_PER_SEC));
#else
	printf_s("Use Time: %f \n", (dur / CLOCKS_PER_SEC));
#endif

	getchar();
	return 1;
}