PS_Method_CCS_test/Pseudo_Spectral_Method/GPM.h

#ifndef GPM_H
#define GPM_H
//代码
//////////////////////*头文件包含部分*/////////////////////////////////////
#include "program.hpp"
#include <cassert>
#include <iostream>
#include <math.h>
#include <armadillo>
using namespace arma;
using namespace std;
using namespace Ipopt;
////////////////////////*类声明部分*////////////////////////////////////////

//////*input结构体部分*//////

/* 过程中的梯度信息结构体 */
struct Grad_struct_info0 {
public:
	/*成员变量*/
	arma::rowvec dJdI;
	int dJdend;
	arma::mat dIdz;
	arma::mat dedS;
	arma::mat dedt;
	arma::mat dedp;

};

/* 过程中的雅可比结构信息 */
struct Jacb_struct_info0 {
public:
	/*成员变量*/
	arma::mat dD_dS_L;
	arma::mat dD_dS_R;
	arma::mat dD_dU;
	arma::urowvec dD_dt;
	arma::umat dD_dp;
	arma::mat dP_dS;
	arma::mat dP_dU;
	arma::mat dP_dt;
	arma::umat dP_dp;
	arma::mat dI_dS;
	arma::mat dI_dU;
	arma::mat dI_dt;
	arma::sp_mat dI_dp;
	arma::mat dE_dS;
	arma::mat dE_dU;
	arma::mat dE_dt;
	arma::mat dE_dp;
};


/* struct结构体 input 部分 */
struct phase0 {
public:
	/*成员变量*/
	arma::vec time;
	arma::mat state;
	arma::mat control;
	arma::rowvec integral;
	arma::vec interval_points;
	arma::mat path;
};

struct guess0 {
public:
	/*成员变量*/
	struct phase0 phase;
};

struct ipoptoptions1 {
public:
	/*成员变量*/
	int print_level;
	char* hessian_approximation;
	char* mu_strategy;
	char* linear_solver;
	double tolerance;
	int maxiterration;
};

struct NLP1 {
public:
	/*成员变量*/
	char *solver;
	struct ipoptoptions1 ipoptoptions;
};

struct phase1 {
public:
	/*成员变量*/
	int inter_num;
	arma::vec colpoints;
	arma::vec i_proportion;
};

struct refinement1 {
public:
	/*成员变量*/
	double Conservative_coefficient;
};

struct mesh1 {
public:
	/*成员变量*/
	char* method;
	double tolerance;
	int maxiterration;
	int colpointsmin;
	int colpointsmax;
	double splitmult;
	double curveratio;
	struct phase1 phase;
	struct refinement1 refinement;
};

struct derivatives1 {
public:
	/*成员变量*/
	char* supplier;
	int derivativelevel;
	char* dependencies;
	int number;
};

struct auxdata1 {
public:
	/*成员变量*/
	double GMe;
	double Sref;
	double M0;
	double C1_v;                                                             
	double Rd_v;                                                            
	double n_v;                                                               
	double m_v;
	double Re;
	double g;
	//新增
	double J2;
	double omige;
	double Ac;
	double ae;
	double re;
	double phif;
	double lamf;
	arma::vec parameter;
};


struct initialtime1 {
public:
	/*成员变量*/
	double lower;
	double upper;
};

struct finaltime1 {
public:
	/*成员变量*/
	double lower;
	double upper;
};

struct state1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct initialstate1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct finalstate1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct control1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct path1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct integral1 {
public:
	/*成员变量*/
	arma::rowvec lower;
	arma::rowvec upper;
};

struct eventgroup1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct event1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct parameter1 {
public:
	/*成员变量*/
	arma::vec lower;
	arma::vec upper;
};

struct phase2 {
public:
	/*成员变量*/
	struct initialtime1 initialtime;
	struct finaltime1 finaltime;
	struct state1 state;
	struct initialstate1 initialstate;
	struct finalstate1 finalstate;
	struct control1 control;
	struct path1 path;
	struct integral1 integral;
	struct eventgroup1 eventgroup;
	struct event1 event;
	struct parameter1 parameter;
};

struct bounds1 {
public:
	/*成员变量*/
	struct phase2 phase;
};

struct interpX{
public:
	/*成员变量*/
	char* method;
};

struct feature1 {
public:
	/*成员变量*/
	int FSE;
};

struct setup0 {
public:
	/*成员变量*/
	int testmode;
	double warmstart;
	struct feature1 feature;
	struct interpX interp;
	char *method;
	struct guess0 guess;
	struct NLP1 NLP;
	struct mesh1 mesh;
	struct derivatives1 derivatives;
	struct auxdata1 auxdata;
	struct bounds1 bounds;
	int control_length;
	int state_length;
	int parameter_length;
	int scales;	
	char* name;
	int Scal;
};
struct phase_result {
public:
	/*成员变量*/
	arma::vec initialstate;
	arma::vec finalstate;
	double initialtime;
	double finaltime;
	arma::rowvec integral;
};
struct input0 {
public:
	/*成员变量*/
	char *name;
	struct setup0 setup;
	struct phase_result phase;
};

/////////////////////////////
struct functions1 {
public:
	/*成员变量*/
	int(*continous1)(mat*, mat*, auxdata1, mat*, mat*, mat*);
	int(*continous2)(cx_mat*, mat*, auxdata1, mat*, mat*, mat*);
	int(*continous3)(mat*, cx_mat*, auxdata1, mat*, mat*, mat*);
	int(*endpoint)(input0*, double*, mat*);
};

////////*output结构体*///////

/* struct结构体 output 部分 */
struct result0 {
public:
	/*成员变量*/
	arma::vec time;
	arma::mat state;
	arma::mat control;
	double objective;
	arma::rowvec mesh_N;
	arma::mat mesh_time;
	auxdata1 auxdata;
	double max_re_error;
};

struct error0 {
public:
	/*成员变量*/
	double max_re_error_max;
	arma::vec re_error_max;
};

struct output0 {
public:
	/*成员变量*/
	struct result0 result;
	arma::vec re_error_max;
	arma::mat Final_State_Estimation;
};

///////////////////////////*函数声明部分*////////////////////////////////////////
int Lagrange_base(mat* point, vec* x, vec* f);
int DLagrange_point(mat* point1, int *k0, int *j, double *f);
double D3Lagrange_point(mat point1, int k0, int j);
int D3L_interval_LGR(mat point1, int Nk, mat* f);
int DL_interval_LGR(mat* point1, int* Nk, mat* f);
int Legendre(int* n, double* x0, double* f);
int Legendre_vector(int* n, vec* x0, vec* f);
int Sum_ZeroPoints(vec* L0, vec* x0, int* Num, mat* Pos);
int Dichotomy_LGR(int* n, vec* Pos, double* tn);
int LGR_points_pro(vec* tn, int* n);
int Weight_LGR_pro(int* N, vec* x, vec* w);
int Intput_Initial(input0* input);
output0* PM_Process(input0* input, functions1 functions, int *iteration_final);
vec Normal(vec b, double bmax, double bmin);
mat Points_Choose(char* method);
int Mesh_initial(int colpointsmax, int interval_num, vec Npoints, mat Col_Points, vec interval_points_ratio,
	mat* MeshGrid, mat* Timestart_Timeend);
int Dismission_check(input0* input, int* Sign_start);
int Method2Nk(const char* method, uvec N, uvec* Nk);
int N_Col2interp(const char* method, uvec N_Col, uvec* N_interp);
int D_pro(mat* MeshGrid, int interval_num, uvec* N_interp, uvec* N_Col, mat* D_phase);
int D_M2Cell(mat* D_phase, int interval_num, uvec* N_interp, field<mat>* D_phasei);
vec Re_Normal(vec b, double bmax, double bmin);
int UY_pro(mat MeshGrid, mat U_old1, mat Y_old1, vec T_old1, mat Ts1_Te1, uvec N_interp, int inter_num, char* method,
	mat* U_phase, mat* Y_phase);
int W_pro(char* method, uvec* N, uvec* N_Col, int interval_num, mat* W_phase);
int N_initial(vec N1, uvec* N, int interval_num);
int Z_pro(mat Y_phase1, mat U_phase1, vec parameter, uvec last_interp, int Col_total, int length_state, int length_control, double tf, double t0, vec* Z);
int Integral_ipopt(vec Z1, mat deta_t_in_col, int Col_total, auxdata1 auxdata,mat W_phase, int length_state, int length_control, vec* integral, int length_parameter, 
					int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path));
int Grad_struct_solve(vec Z0, mat deta_t_in_col, int Col_total, auxdata1 auxdata, mat W_phase, int length_state,
	int length_control, int length_integral, Grad_struct_info0* Grad_struct_info, int length_parameter,
	int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
	int(*Objective_ptr)(input0* input, double* output, mat* event));
int End_Objective_ipopt(vec Z1, int Col_total, int length_state, int length_integral, int length_parameter, int length_control, 
						int(*Objective_ptr)(input0* input, double* output, mat* event),
						double* J_end);
int Istruct(vec Z1, int Col_total, int length_state, int length_control, int length_integral, mat W_phase, auxdata1 auxdata, 
	int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
	vec* I);
int IZstruct(vec Z0, int Col_total, int length_state, int length_control, int length_integral,mat W_phase, auxdata1 auxdata, sp_mat*I_Z_struct0, 
			int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path));
int Fmin_pro(input0 input, int Col_total, int length_state, int length_path, int length_integral, int length_event, vec* Fmin);
int Fmax_pro(input0 input, int Col_total, int length_state, int length_path, int length_integral, int length_event, vec* Fmax);
int Zmin_pro(input0 input, int Col_total, vec* Zmin);
int Zmax_pro(input0 input, int Col_total, vec* Zmax);
int F_pro_struct(vec Z1, mat D_phase, uvec first_col, mat deta_t_in_col, int length_path, int length_parameter, int Col_total, auxdata1 auxdata,
				mat W_phase, int z_len, int length_event, int length_state, int length_control, int length_integral, int interval_num, uvec last_interp,
				int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
				int(*Objective_ptr)(input0* input, double* output, mat* event),
				vec* F1);
int Jacb_struct_pro0(mat Z0, mat D_phase, uvec first_col, mat deta_t_in_col, int length_event, int length_path, int length_parameter,
				int Col_total, auxdata1 auxdata, mat W_phase, int length_state, int length_control, int length_integral, int length_F0,
				int interval_num, uvec last_interp, 
				int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
				int(*Objective_ptr)(input0* input, double* output, mat* event),
				sp_mat *Jacb_struct0);
int Jacb_solve(sp_mat Jacb_struct01, int Col_total, int length_event, int length_path, int length_state, int length_control, int length_integral, int length_parameter, Jacb_struct_info0* Jacb_struct_info);
int vector_ARM2NLP(vec vector_ARM, Number* vector_NLP, int n);
int vector_NLP2ARM(const Number* vector_NLP, vec* vector_ARM, int n);
int Objective_ipopt(vec Z1, mat deta_t_in_col, int Col_total, auxdata1 auxdata, mat W_phase, int length_state, int length_control, int length_parameter,
					int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
					int(*Objective_ptr)(input0* input, double* output, mat* event),
					double* J);
int Grad_pro(vec Z0, Grad_struct_info0 Grad_struct_info, mat deta_t_in_col, int Col_total, auxdata1 auxdata, mat W_phase, int length_state, 
			int length_control, int length_integral, int length_parameter, double eps1, char *deri_supplier,
			int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
			int(*Objective_ptr)(input0* input, double* output, mat* event),
			sp_mat *Grad);
int F_pro(vec Z1, mat D_phase, uvec first_col, mat deta_t_in_col, int length_path, int Col_total, auxdata1 auxdata, mat W_phase, 
		int z_len, int length_event, int length_state, int length_control, int length_integral, int length_parameter, 
		int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
		int(*Objective_ptr)(input0* input, double* output, mat* event),
		vec* F1);
int Jacb_pro(vec Z0, field<mat> D_phasei, uvec N_Col, mat deta_t_in_col, sp_mat I_Z_struct0, Jacb_struct_info0 Jacb_struct_info, int length_F0, int length_path, int Col_total, auxdata1 auxdata,
	mat W_phase, int length_event, int length_state, int length_control, int length_integral, int length_parameter, int interval_num, double eps1, char *deri_supplier, int length_Z0, 
	int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
	int(*Objective_ptr)(input0* input, double* output, mat* event),
	sp_mat* Jacb);
int sp_ARM2NLP(sp_mat X, Index* iRow, Index *jCol, int nnz_jac);
int values_initial(Number* values, Index nnz_jac);
int values_set(sp_mat Jacb, sp_mat Jacb_struct0, Number* values, Index nnz_jac);
void get_var_num(int varnum0);
void get_constraint_num(int constraint_num0);
void get_nnz_jac(int nnz_jac0);
void get_nnz_h(int nnz_h0);
void get_x0(vec x00);
void get_bounds(vec Fmin0, vec Fmax0, vec Zmin0, vec Zmax0);
void get_data(auxdata1 auxdata00, mat W_phase00, mat deta_t_in_col00, int Col_total00, int length_state00, int length_control00,
	Grad_struct_info0 Grad_struct_info00, int length_integral00, double eps00, mat D_phase00, uvec first_col00, int length_path00,
	int length_event00, sp_mat Jacb_struct00, field<mat> D_phasei0, uvec N_Col00, sp_mat I_Z_struct00, int interval_num00,
	int length_Z00, int length_parameter00, Jacb_struct_info0 Jacb_struct_info00, char *deri_supplier00, int length_F000, 
	int(*Dynamics_ptr00)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
	int(*Objective_ptr)(input0* input, double* output, mat* event));
int ipopt_apply(int var_num, int constraint_num, int nnz_jac, int nnz_h, vec x00, vec Fmin0, vec Fmax0, vec Zmin0, vec Zmax0,
	auxdata1 auxdata00, mat W_phase00, mat deta_t_in_col00, int Col_total00, int length_state00, int length_control00,
	Grad_struct_info0 Grad_struct_info00, int length_integral00, double eps00, mat D_phase00, uvec first_col00, int length_path00,
	int length_event00, sp_mat Jacb_struct00, field<mat> D_phasei0, uvec N_Col00, sp_mat I_Z_struct00, int interval_num00,
	int length_Z00, int length_parameter00, Jacb_struct_info0 Jacb_struct_info00, char *deri_supplier00, int length_F000, NLP1 NLP,
	int(*Dynamics_ptr00)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
	int(*Objective_ptr)(input0* input, double* output, mat* event),
	vec* z_ans, vec* zl_temp, vec* zu_temp, vec* lambda_temp, int* iter_num, double *cpu_time);
int Status_Display(int status);
int Z2YU_interval(vec x_ans, int Col_total, int length_state, int length_control, mat *Y_phase_ans, mat *U_phase_ans);
mat Re_Normal(mat b, double bmax, double bmin);
int Time_ID(vec* Ts1_Te1, double t0, double tf, double t0_ans, double tf_ans);
int Time2Newtime(mat Ts1_Te1, mat MeshGrid, double max_time, double min_time, double t0_ans, double tf_ans, int interval_num, int Col_total, uvec last_col, vec* Time_interval_ans);
mat Normal(mat b, double bmax, double bmin);
int Jacb_struct_Re(Jacb_struct_info0 Jacb_struct_info, mat MeshGrid_nextstep, int length_state, int length_control, int length_path, int length_event, int length_integral, int length_parameter, sp_mat* Jacb_struct1);
// 终端误差状态估计
int FinalStateEstimate(mat U_ans, mat Y_ans, mat MeshGrid, double t0_ans, double tf_ans, uvec N, char *method, vec time_ans,
						int interval_num, mat Ts1_Te1_ans, auxdata1 auxdata, int length_state, int length_control, char *interp_method,
						int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path),
						mat* FSE);
// 结果评价与网格细化策略
int result_assess(arma::mat* MeshGrid_nextstep, arma::mat* Ts1_Te1_nextstep, error0* error, arma::mat* inum_interval_increase,
	arma::mat U_ans, arma::mat Y_ans, arma::mat MeshGrid, arma::mat time1_ans, double t0_ans, double tf_ans,
	arma::uvec N, mesh1 mesh, arma::mat Col_Points, const char* method,
	unsigned int interval_num, arma::mat Ts1_Te1_ans, auxdata1 auxdata,
	unsigned int length_state, unsigned length_control, int niu, int print_level, int testmode, double CC, 
	int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path));

mat Mesh_choose(mesh1 mesh, arma::mat Col_Points, arma::uvec N);		// 根据阶数进行网格信息的选择
arma::field<arma::mat> D_pro1(arma::mat MeshGrid, arma::uword interval_num);// 更新阶段内的左动力学微分矩阵D
arma::field<arma::mat> U_pro(arma::mat MeshGridnew, arma::mat MeshGrid, arma::mat U_old, arma::uvec Nk, arma::uword inter_num, unsigned length_control, int deta_N); // 更新控制量矩阵U
arma::field<arma::mat> Y_pro(arma::mat MeshGridnew, arma::mat MeshGrid, arma::mat Y_old, arma::uvec Nk, arma::uword inter_num, unsigned length_state, int deta_N); // 更新状态量矩阵Y
arma::field<arma::mat> Integral_Matrix(arma::field<arma::mat> D_phase_ans, arma::uword interval_num);	// 计算拉格朗日隐式积分矩阵
field<mat> First_deriv_approximation(field<mat> U_phase_ans, field<mat> Y_phase_ans, uword interval_num, auxdata1 auxdata, 
									int(*Dynamics_ptr)(mat* StateVar, mat* ControlVar, auxdata1 auxdata, mat* dState, mat* integrand, mat* path));


// 新增函数部分
int Get_Rou_Mach(vec H, vec V, vec* Rou, vec* Ma);
int alpha2LD(vec alpha_rad, vec Rou, vec Ma, vec V, auxdata1 auxdata, vec* L, vec* D);

#endif
#pragma once