1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397339833993400340134023403340434053406340734083409341034113412341334143415341634173418341934203421342234233424342534263427342834293430343134323433343434353436343734383439344034413442344334443445344634473448344934503451345234533454345534563457345834593460346134623463346434653466346734683469347034713472347334743475347634773478347934803481348234833484348534863487348834893490349134923493349434953496349734983499350035013502350335043505350635073508350935103511351235133514351535163517351835193520352135223523352435253526352735283529353035313532353335343535353635373538353935403541354235433544354535463547354835493550355135523553355435553556355735583559356035613562356335643565356635673568356935703571357235733574357535763577357835793580358135823583358435853586358735883589359035913592359335943595359635973598359936003601360236033604360536063607360836093610361136123613361436153616361736183619362036213622362336243625362636273628362936303631363236333634363536363637363836393640 |
- #include "cmscp_socpInfo.h"
- #include "freeSocpInfo.h"
- void socpInfoSparseDim(cmscp_socpInfo** socpInfo0, cmscp_consTypeInfo* consTypeInfo, cmscp_sparseInfo* sparseInfo, cmscp_setup* setup)
- {
- // 函数功能:获取二阶锥规划问题的参数的维数以及分配相应的空间
- // 输入:
- // socpInfo0:无维数和未分配二阶锥规划问题参数信息
- // consTypeInfo:约束类型信息
- // sparseInfo:稀疏(非零元位置)信息
- // setup:问题设置参数
- // 输出:
- // socpInfo0:有维数和已分配二阶锥规划问题参数信息
- // 假设函数的自变量按照[x, u, p]排序
- // 假设优化变量的排列顺序为
- // [{x1}, { u1 }, t01, tf1, p1, Tr1, .., { xi }, { ui }, t0i, tfi, pi, Tri, ..., { xN }, { uN }, t0N, tfN, pN, TrN];
- // 其中{ xi } = { xi_1,...,xi_j,...,xi_n }, xi_j表示第i段的第j个状态向量
- // { ui } = { ui_1,...,ui_j,...,ui_n }, xi_j表示第i段的第j个控制输入
- // t0i, tfi表示第i段初始和终端时刻
- // pi表示第i段参数变量
- // Tri表示软信赖域边界
- // 谢磊:2022 / 6 / 27编写
- // 变量声明
- int numPhase, numOptVar, * numPhaseOptVar, iphase, numMesh, numState, numControl, numParameter, numRowAEqua, numNzAEqua, numRowAIneq, numNzAIneq, numRowASoc, numNzASoc, numSubASoc, numNlPath, numSubSocPath, * dimSocPath, * contStateNzPtrRow, * contControlNzPtrRow, * contParameterNzPtrRow, numNzPhaseAEqua, numRowPhaseAIneq, numRowPhaseAEqua, numNzPhaseAIneq, iType, numType, * idxType, cnt, idxRow, numPhaseSubASoc, numRowPointASoc, numRowPhaseASoc, idxStartRowSoc, idxEndRowSoc, numNzPhaseASoc, numNlEndp, numSubSocEndp, * dimSocEndp, * endpIniStateNzPtrRow, * endpFinStateNzPtrRow, * endpIniTimeNzPtrRow, * endpFinTimeNzPtrRow, * endpParameterNzPtrRow, numNlLink, numSubSocLink, * dimSocLink, * linkLeftStatefNzPtrRow, * linkLeftTfNzPtrRow, * linkLeftParameterNzPtrRow, * linkRightState0NzPtrRow, * linkRightT0NzPtrRow, * linkRightParameterNzPtrRow, * objEventIniStateNzPos, * objEventFinStateNzPos, * objEventIniTimeNzPos, * objEventFinTimeNzPos, * objEventParameterNzPos;
- double* c;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_nzPos_cont** contNzPos;
- cmscp_nzPos_endp** objEventNzPos;
- cmscp_nzPos_endp** endpNzPos;
- cmscp_nzPos_link** linkNzPos;
- cmscp_consType** endpConsType;
- cmscp_consType** pathConsType;
- cmscp_consType_bound* boundConsType;
- cmscp_consType** linkConsType;
- cmscp_consType* eventConsType;
- cmscp_consType** boundIniTimeType;
- cmscp_consType** boundFinTimeType;
- cmscp_consType** boundIniStateType;
- cmscp_consType** boundFinStateType;
- cmscp_consType** boundStateType;
- cmscp_consType** boundControlType;
- cmscp_consType** boundParameterType;
- cmscp_socpInfo* socpInfo;
- // 1. 基本参数
- dimension = setup->dimension; // 问题维数信息
- mesh = setup->mesh; // 网格信息
- numPhase = setup->numPhase; // 阶段数
- // 2. 优化变量的个数
- numOptVar = 0;
- numPhaseOptVar = (int*)malloc(numPhase * sizeof(int));
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- numPhaseOptVar[iphase] = numMesh * (numState + numControl) + numParameter + 3; // 保留软信赖域边
- numOptVar = numOptVar + numPhaseOptVar[iphase];
- }
- // 3. 目标函数系数
- c = (double*)malloc(numOptVar * sizeof(double));
- // 4. 统计三类矩阵的非零元个数和行数:不等式约束矩阵AIneq、等式约束矩阵AEqua和二阶锥约束矩阵ASoc
- numRowAEqua = 0;// 等式约束矩阵的行数
- numNzAEqua = 0;// 等式约束矩阵的非零个数
- numRowAIneq = 0;// 不等式约束矩阵的行数
- numNzAIneq = 0;// 不等式约束矩阵的非零个数
- numRowASoc = 0;// 二阶锥矩阵的行数
- numNzASoc = 0; // 二阶锥约束矩阵的非零元
- numSubASoc = 0;// 二阶锥矩阵中子二阶锥的个数
- contNzPos = sparseInfo->contNzPos;// 连续函数雅可比稀疏信息
- endpNzPos = sparseInfo->endpNzPos;// 端点函数雅可比稀疏信息
- objEventNzPos = sparseInfo->objEventNzPos;// 目标函数雅可比稀疏信息
- linkNzPos = sparseInfo->linkNzPos;// 衔接函数雅可比稀疏信息
- // 约束类型
- endpConsType = consTypeInfo->endp;
- linkConsType = consTypeInfo->link;
- pathConsType = consTypeInfo->path;
- boundConsType = consTypeInfo->bound;
- boundIniTimeType = boundConsType->iniTime;
- boundFinTimeType = boundConsType->finTime;
- boundIniStateType = boundConsType->iniState;
- boundFinStateType = boundConsType->finState;
- boundStateType = boundConsType->state;
- boundControlType = boundConsType->control;
- boundParameterType = boundConsType->parameter;
- eventConsType = consTypeInfo->event;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 维数信息
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 连续函数部分
- // 约束维度信息
- numNlPath = dimension->phase[iphase]->path->nl;
- numSubSocPath = dimension->phase[iphase]->path->soc->n;
- dimSocPath = dimension->phase[iphase]->path->soc->q;
- // 连续函数雅可比矩阵稀疏信息
- contStateNzPtrRow = contNzPos[iphase]->state->ptrRow;
- contControlNzPtrRow = contNzPos[iphase]->control->ptrRow;
- contParameterNzPtrRow = contNzPos[iphase]->parameter->ptrRow;
- // 动力学方程部分
- numNzPhaseAEqua = 0;
- numNzPhaseAEqua = numNzPhaseAEqua + 2 * contStateNzPtrRow[numState];// 状态变量部分
- numNzPhaseAEqua = numNzPhaseAEqua + 2 * contControlNzPtrRow[numState];// 控制变量部分
- numNzPhaseAEqua = numNzPhaseAEqua + contParameterNzPtrRow[numState];// 静态参数部分
- numNzPhaseAEqua = numNzPhaseAEqua + numState * 2;// 初始和终端时刻部分
- numNzPhaseAEqua = (numMesh - 1) * numNzPhaseAEqua;
- numNzAEqua = numNzAEqua + numNzPhaseAEqua;
- numRowAEqua = numRowAEqua + (numMesh - 1) * numState;
- // 路径约束部分
- // 不等式和等式部分
- numRowPhaseAIneq = 0;// 各段不等式约束的行数
- numRowPhaseAEqua = 0;// 各段等式约束的行数
- numNzPhaseAIneq = 0;// 各段不等式约束的非零元
- numNzPhaseAEqua = 0;// 各段等式约束的非零元
- for (iType = 0; iType < 4; iType++)
- {
- numType = pathConsType[iphase]->type[iType].num;
- idxType = pathConsType[iphase]->type[iType].idx;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt] + numState; // 在连续函数雅可比矩阵中行位置
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (contStateNzPtrRow[idxRow + 1] - contStateNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (contControlNzPtrRow[idxRow + 1] - contControlNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (contParameterNzPtrRow[idxRow + 1] - contParameterNzPtrRow[idxRow]);
- }
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt] + numState;
- numNzPhaseAEqua = numNzPhaseAEqua + contStateNzPtrRow[idxRow + 1] - contStateNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + contControlNzPtrRow[idxRow + 1] - contControlNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + contParameterNzPtrRow[idxRow + 1] - contParameterNzPtrRow[idxRow];
- }
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt] + numState;
- numNzPhaseAIneq = numNzPhaseAIneq + contStateNzPtrRow[idxRow + 1] - contStateNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + contControlNzPtrRow[idxRow + 1] - contControlNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + contParameterNzPtrRow[idxRow + 1] - contParameterNzPtrRow[idxRow];
- }
- }
- }
- }
- numRowPhaseAIneq = numMesh * numRowPhaseAIneq;
- numRowPhaseAEqua = numMesh * numRowPhaseAEqua;
- numRowAIneq = numRowAIneq + numRowPhaseAIneq;
- numRowAEqua = numRowAEqua + numRowPhaseAEqua;
- numNzPhaseAIneq = numMesh * numNzPhaseAIneq;
- numNzPhaseAEqua = numMesh * numNzPhaseAEqua;
- numNzAIneq = numNzAIneq + numNzPhaseAIneq;
- numNzAEqua = numNzAEqua + numNzPhaseAEqua;
- if (numSubSocPath)
- {
- // 二阶锥个数
- numPhaseSubASoc = numMesh * numSubSocPath;// 各段二阶锥矩阵个数
- numSubASoc = numSubASoc + numPhaseSubASoc;// 总二阶锥个数
- // 二阶锥矩阵行数
- numRowPointASoc = cmscp_sum(dimSocPath, numSubSocPath);// 单个离散点处的二阶锥矩阵行数
- numRowPhaseASoc = numMesh * numRowPointASoc;// 各段二阶锥矩阵行数
- numRowASoc = numRowASoc + numMesh * numRowPhaseASoc;// 总二阶锥矩阵行数
- // 单段二阶锥矩阵的非零元个数
- idxStartRowSoc = numState + numNlPath;// 二阶锥矩阵在连续函数雅可比矩阵中的起始行
- idxEndRowSoc = numState + numNlPath + numRowPointASoc;// 二阶锥矩阵在连续函数雅可比矩阵中的结束行
- numNzPhaseASoc = numNzPhaseASoc + contStateNzPtrRow[idxEndRowSoc + 1] - contStateNzPtrRow[idxStartRowSoc + 1];
- numNzPhaseASoc = numNzPhaseASoc + contControlNzPtrRow[idxEndRowSoc + 1] - contControlNzPtrRow[idxStartRowSoc + 1];
- numNzPhaseASoc = numNzPhaseASoc + contParameterNzPtrRow[idxEndRowSoc + 1] - contParameterNzPtrRow[idxStartRowSoc + 1];
- numNzPhaseASoc = numMesh * numNzPhaseASoc;
- // 总二阶锥矩阵非零元个数
- numNzASoc = numNzASoc + numNzPhaseASoc;
- }
- //// 端点函数部分
- // 约束维度信息
- numNlEndp = dimension->phase[iphase]->endpoint->nl;
- numSubSocEndp = dimension->phase[iphase]->endpoint->soc->n;
- dimSocEndp = dimension->phase[iphase]->endpoint->soc->q;
- // 端点函数雅可比矩阵稀疏信息
- endpIniStateNzPtrRow = endpNzPos[iphase]->iniState->ptrRow;
- endpFinStateNzPtrRow = endpNzPos[iphase]->finState->ptrRow;
- endpIniTimeNzPtrRow = endpNzPos[iphase]->iniTime->ptrRow;
- endpFinTimeNzPtrRow = endpNzPos[iphase]->finTime->ptrRow;
- endpParameterNzPtrRow = endpNzPos[iphase]->parameter->ptrRow;
- // 不等式和等式部分
- numRowPhaseAIneq = 0;// 各段不等式约束的行数
- numRowPhaseAEqua = 0;// 各段等式约束的行数
- numNzPhaseAIneq = 0;// 各段不等式约束的非零元
- numNzPhaseAEqua = 0;// 各段等式约束的非零元
- for (iType = 0; iType < 4; iType++)
- {
- numType = endpConsType[iphase]->type[iType].num;
- idxType = endpConsType[iphase]->type[iType].idx;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (endpIniStateNzPtrRow[idxRow + 1] - endpIniStateNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (endpFinStateNzPtrRow[idxRow + 1] - endpFinStateNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (endpIniTimeNzPtrRow[idxRow + 1] - endpIniTimeNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (endpFinTimeNzPtrRow[idxRow + 1] - endpFinTimeNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (endpParameterNzPtrRow[idxRow + 1] - endpParameterNzPtrRow[idxRow]);
- }
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAEqua = numNzPhaseAEqua + endpIniStateNzPtrRow[idxRow + 1] - endpIniStateNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + endpFinStateNzPtrRow[idxRow + 1] - endpFinStateNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + endpIniTimeNzPtrRow[idxRow + 1] - endpIniTimeNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + endpFinTimeNzPtrRow[idxRow + 1] - endpFinTimeNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + endpParameterNzPtrRow[idxRow + 1] - endpParameterNzPtrRow[idxRow];
- }
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAIneq = numNzPhaseAIneq + endpIniStateNzPtrRow[idxRow + 1] - endpIniStateNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + endpFinStateNzPtrRow[idxRow + 1] - endpFinStateNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + endpIniTimeNzPtrRow[idxRow + 1] - endpIniTimeNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + endpFinTimeNzPtrRow[idxRow + 1] - endpFinTimeNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + endpParameterNzPtrRow[idxRow + 1] - endpParameterNzPtrRow[idxRow];
- }
- }
- }
- }
- numRowAIneq = numRowAIneq + numRowPhaseAIneq;
- numRowAEqua = numRowAEqua + numRowPhaseAEqua;
- numNzAIneq = numNzAIneq + numNzPhaseAIneq;
- numNzAEqua = numNzAEqua + numNzPhaseAEqua;
- // 二阶锥约束部分
- if (numSubSocEndp)
- {
- // 二阶锥个数
- numSubASoc = numSubASoc + numSubSocEndp; // 总二阶锥个数, 每段就一个端点约束函数
- // 二阶锥矩阵行数
- numRowPhaseASoc = cmscp_sum(dimSocEndp, numSubSocEndp); // 各段二阶锥矩阵行数
- numRowASoc = numRowASoc + numRowPhaseASoc; // 总二阶锥矩阵行数
- // 单段矩阵的非零元个数
- idxStartRowSoc = numNlEndp; // 二阶锥矩阵在连续函数雅可比矩阵中的起始行
- idxEndRowSoc = numNlEndp + numRowPhaseASoc; // 二阶锥矩阵在连续函数雅可比矩阵中的结束行
- numNzPhaseASoc = 0;
- numNzPhaseASoc = numNzPhaseASoc + endpIniStateNzPtrRow[idxEndRowSoc] - endpIniStateNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + endpFinStateNzPtrRow[idxEndRowSoc] - endpFinStateNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + endpIniTimeNzPtrRow[idxEndRowSoc] - endpIniTimeNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + endpFinTimeNzPtrRow[idxEndRowSoc] - endpFinTimeNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + endpParameterNzPtrRow[idxEndRowSoc] - endpParameterNzPtrRow[idxStartRowSoc];
- // 总二阶锥矩阵非零元个数
- numNzASoc = numNzASoc + numNzPhaseASoc;
- }
- if (iphase < numPhase - 1)
- {
- // 约束维数信息
- numNlLink = dimension->phase[iphase]->linkage->nl;// 非线性部分维数
- numSubSocLink = dimension->phase[iphase]->linkage->soc->n;
- dimSocLink = dimension->phase[iphase]->linkage->soc->q;
- // 稀疏信息
- linkLeftStatefNzPtrRow = linkNzPos[iphase]->leftStatef->ptrRow;
- linkLeftTfNzPtrRow = linkNzPos[iphase]->leftTf->ptrRow;
- linkLeftParameterNzPtrRow = linkNzPos[iphase]->leftParameter->ptrRow;
- linkRightState0NzPtrRow = linkNzPos[iphase]->rightState0->ptrRow;
- linkRightT0NzPtrRow = linkNzPos[iphase]->rightT0->ptrRow;
- linkRightParameterNzPtrRow = linkNzPos[iphase]->rightParameter->ptrRow;
- // 统计各类非线性约束矩阵的行数
- numRowPhaseAIneq = 0;// 各段不等式约束的行数
- numRowPhaseAEqua = 0;// 各段等式约束的行数
- numNzPhaseAIneq = 0;// 各段不等式约束的非零元
- numNzPhaseAEqua = 0;// 各段等式约束的非零元
- for (iType = 0; iType < 4; iType++)
- {
- numType = linkConsType[iphase]->type[iType].num;
- idxType = linkConsType[iphase]->type[iType].idx;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkLeftStatefNzPtrRow[idxRow + 1] - linkLeftStatefNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkLeftTfNzPtrRow[idxRow + 1] - linkLeftTfNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkLeftParameterNzPtrRow[idxRow + 1] - linkLeftParameterNzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkRightState0NzPtrRow[idxRow + 1] - linkRightState0NzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkRightT0NzPtrRow[idxRow + 1] - linkRightT0NzPtrRow[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (linkRightParameterNzPtrRow[idxRow + 1] - linkRightParameterNzPtrRow[idxRow]);
- }
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAEqua = numNzPhaseAEqua + linkLeftStatefNzPtrRow[idxRow + 1] - linkLeftStatefNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + linkLeftTfNzPtrRow[idxRow + 1] - linkLeftTfNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + linkLeftParameterNzPtrRow[idxRow + 1] - linkLeftParameterNzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + linkRightState0NzPtrRow[idxRow + 1] - linkRightState0NzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + linkRightT0NzPtrRow[idxRow + 1] - linkRightT0NzPtrRow[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + linkRightParameterNzPtrRow[idxRow + 1] - linkRightParameterNzPtrRow[idxRow];
- }
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- for (cnt = 0; cnt < numType; cnt++)
- {
- idxRow = idxType[cnt];
- numNzPhaseAIneq = numNzPhaseAIneq + linkLeftStatefNzPtrRow[idxRow + 1] - linkLeftStatefNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + linkLeftTfNzPtrRow[idxRow + 1] - linkLeftTfNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + linkLeftParameterNzPtrRow[idxRow + 1] - linkLeftParameterNzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + linkRightState0NzPtrRow[idxRow + 1] - linkRightState0NzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + linkRightT0NzPtrRow[idxRow + 1] - linkRightT0NzPtrRow[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + linkRightParameterNzPtrRow[idxRow + 1] - linkRightParameterNzPtrRow[idxRow];
- }
- }
- }
- }
- numRowAIneq = numRowAIneq + numRowPhaseAIneq;
- numRowAEqua = numRowAEqua + numRowPhaseAEqua;
- numNzAIneq = numNzAIneq + numNzPhaseAIneq;
- numNzAEqua = numNzAEqua + numNzPhaseAEqua;
- if (numSubSocLink)
- {
- // 二阶锥个数
- numSubASoc = numSubASoc + numSubSocLink;// 总二阶锥个数, 每段就一个端点约束函数
- // 二阶锥矩阵行数
- numRowPhaseASoc = cmscp_sum(dimSocLink, numSubSocLink);// 各段二阶锥矩阵行数
- numRowASoc = numRowASoc + numRowPhaseASoc;// 总二阶锥矩阵行数
- // 单段矩阵的非零元个数
- idxStartRowSoc = numNlLink;// 二阶锥矩阵在连续函数雅可比矩阵中的起始行
- idxEndRowSoc = numNlLink + numRowPhaseASoc;// 二阶锥矩阵在连续函数雅可比矩阵中的结束行
- numNzPhaseASoc = 0;
- numNzPhaseASoc = numNzPhaseASoc + linkLeftStatefNzPtrRow[idxEndRowSoc] - linkLeftStatefNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + linkLeftTfNzPtrRow[idxEndRowSoc] - linkLeftTfNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + linkLeftParameterNzPtrRow[idxEndRowSoc] - linkLeftParameterNzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + linkRightState0NzPtrRow[idxEndRowSoc] - linkRightState0NzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + linkRightT0NzPtrRow[idxEndRowSoc] - linkRightT0NzPtrRow[idxStartRowSoc];
- numNzPhaseASoc = numNzPhaseASoc + linkRightParameterNzPtrRow[idxEndRowSoc] - linkRightParameterNzPtrRow[idxStartRowSoc];
- // 总二阶锥矩阵非零元个数
- numNzASoc = numNzASoc + numNzPhaseASoc;
- }
- }
- //// 上下边界约束
- // 各段不等式和等式约束的行数也等于非零元个数
- numRowPhaseAEqua = 0;
- numRowPhaseAIneq = 0;
- // 统计行数
- for (iType = 0; iType < 4; iType++)
- {
- // 初始时刻变量
- numType = boundIniTimeType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- }
- }
- // 终端时刻变量
- numType = boundFinTimeType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- }
- }
- // 初始状态变量
- numType = boundIniStateType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- }
- }
- // 终端状态变量
- numType = boundFinStateType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- }
- }
- // 过程状态变量
- numType = boundStateType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * (numMesh - 2) * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + (numMesh - 2) * numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + (numMesh - 2) * numType;
- }
- }
- // 控制变量
- numType = boundControlType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numMesh * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numMesh * numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numMesh * numType;
- }
- }
- // 静态参数变量
- numType = boundParameterType[iphase]->type[iType].num;
- if (numType)
- {
- if (iType == 0)
- {
- numRowPhaseAIneq = numRowPhaseAIneq + 2 * numType;
- }
- else if (iType == 1)
- {
- numRowPhaseAEqua = numRowPhaseAEqua + numType;
- }
- else
- {
- numRowPhaseAIneq = numRowPhaseAIneq + numType;
- }
- }
- }
- numRowAIneq = numRowAIneq + numRowPhaseAIneq;
- numRowAEqua = numRowAEqua + numRowPhaseAEqua;
- numNzAIneq = numNzAIneq + numRowPhaseAIneq;
- numNzAEqua = numNzAEqua + numRowPhaseAEqua;
- //// 信赖域约束
- // 二阶锥约束个数
- // 状态:numMesh、控制:numMesh、初始时刻:1、终端时刻:1、静态参数:1
- numSubASoc = numSubASoc + numMesh * 2 + 3;
- // 二阶锥约束矩阵行数和非零元个数相等
- // 每个锥的维数比对应变量多1
- numRowPhaseASoc = numMesh * (numState + numControl + 2) + numParameter + 5;
- numRowASoc = numRowASoc + numRowPhaseASoc;
- numNzASoc = numNzASoc + numRowPhaseASoc;
- }
- // 事件约束
- for (iType = 0; iType < 4; iType++)
- {
- numType = eventConsType->type[iType].num;
- idxType = eventConsType->type[iType].idx;
- if (numType)
- {
- if (iType == 0)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- // 统计非零元素个数
- // 第idxType(cnt)个event约束被分为numPhase段
- idxRow = idxType[cnt] + 1; // 加偏置1,因为稀疏信息是objEvent的,第一行是目标函数obj的
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 稀疏信息
- objEventIniStateNzPos = objEventNzPos[iphase]->iniState->ptrRow;
- objEventFinStateNzPos = objEventNzPos[iphase]->finState->ptrRow;
- objEventIniTimeNzPos = objEventNzPos[iphase]->iniTime->ptrRow;
- objEventFinTimeNzPos = objEventNzPos[iphase]->finTime->ptrRow;
- objEventParameterNzPos = objEventNzPos[iphase]->parameter->ptrRow;
- // 每段非零元个数
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (objEventIniStateNzPos[idxRow + 1] - objEventIniStateNzPos[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (objEventFinStateNzPos[idxRow + 1] - objEventFinStateNzPos[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (objEventIniTimeNzPos[idxRow + 1] - objEventIniTimeNzPos[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (objEventFinTimeNzPos[idxRow + 1] - objEventFinTimeNzPos[idxRow]);
- numNzPhaseAIneq = numNzPhaseAIneq + 2 * (objEventParameterNzPos[idxRow + 1] - objEventParameterNzPos[idxRow]);
- }
- numRowAIneq = numRowAIneq + 1;
- }
- }
- else if (iType == 1)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- // 统计非零元素个数
- // 第idxType(cnt)个event约束被分为numPhase段
- idxRow = idxType[cnt] + 1; // 加偏置1,因为稀疏信息是objEvent的,第一行是目标函数obj的
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 稀疏信息
- objEventIniStateNzPos = objEventNzPos[iphase]->iniState->ptrRow;
- objEventFinStateNzPos = objEventNzPos[iphase]->finState->ptrRow;
- objEventIniTimeNzPos = objEventNzPos[iphase]->iniTime->ptrRow;
- objEventFinTimeNzPos = objEventNzPos[iphase]->finTime->ptrRow;
- objEventParameterNzPos = objEventNzPos[iphase]->parameter->ptrRow;
- // 每段非零元个数
- numNzPhaseAEqua = numNzPhaseAEqua + objEventIniStateNzPos[idxRow + 1] - objEventIniStateNzPos[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + objEventFinStateNzPos[idxRow + 1] - objEventFinStateNzPos[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + objEventIniTimeNzPos[idxRow + 1] - objEventIniTimeNzPos[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + objEventFinTimeNzPos[idxRow + 1] - objEventFinTimeNzPos[idxRow];
- numNzPhaseAEqua = numNzPhaseAEqua + objEventParameterNzPos[idxRow + 1] - objEventParameterNzPos[idxRow];
- }
- numRowAEqua = numRowAEqua + 1;
- }
- }
- else
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- // 统计非零元素个数
- // 第idxType(cnt)个event约束被分为numPhase段
- idxRow = idxType[cnt] + 1; // 加偏置1,因为稀疏信息是objEvent的,第一行是目标函数obj的
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 稀疏信息
- objEventIniStateNzPos = objEventNzPos[iphase]->iniState->ptrRow;
- objEventFinStateNzPos = objEventNzPos[iphase]->finState->ptrRow;
- objEventIniTimeNzPos = objEventNzPos[iphase]->iniTime->ptrRow;
- objEventFinTimeNzPos = objEventNzPos[iphase]->finTime->ptrRow;
- objEventParameterNzPos = objEventNzPos[iphase]->parameter->ptrRow;
- // 每段非零元个数
- numNzPhaseAIneq = numNzPhaseAIneq + objEventIniStateNzPos[idxRow + 1] - objEventIniStateNzPos[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + objEventFinStateNzPos[idxRow + 1] - objEventFinStateNzPos[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + objEventIniTimeNzPos[idxRow + 1] - objEventIniTimeNzPos[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + objEventFinTimeNzPos[idxRow + 1] - objEventFinTimeNzPos[idxRow];
- numNzPhaseAIneq = numNzPhaseAIneq + objEventParameterNzPos[idxRow + 1] - objEventParameterNzPos[idxRow];
- }
- numRowAIneq = numRowAIneq + 1;
- }
- }
- }
- }
- // 分配约束相关空间
- socpInfo = (cmscp_socpInfo*)malloc(sizeof(cmscp_socpInfo));
- socpInfo->bEqua = (double*)malloc(numRowAEqua * sizeof(double));
- socpInfo->bCone = (double*)malloc((numRowAIneq + numRowASoc) * sizeof(double));
- socpInfo->AEqua = (cmscp_crsSpmat*)malloc(sizeof(cmscp_crsSpmat));
- socpInfo->ACone = (cmscp_crsSpmat*)malloc(sizeof(cmscp_crsSpmat));
- socpInfo->AEqua->ptrRow = (int*)malloc((numRowAEqua + 1) * sizeof(int));
- socpInfo->AEqua->ptrRow[numRowAEqua] = numNzAEqua;
- socpInfo->AEqua->idxCol = (int*)malloc(numNzAEqua * sizeof(int));
- socpInfo->AEqua->v = (double*)malloc(numNzAEqua * sizeof(double));
- socpInfo->AEqua->m = numRowAEqua;
- socpInfo->AEqua->n = numOptVar;
- socpInfo->ACone->ptrRow = (int*)malloc((numRowAIneq + numRowASoc + 1) * sizeof(int));
- socpInfo->ACone->ptrRow[numRowAIneq + numRowASoc] = numNzAIneq + numNzASoc;
- socpInfo->ACone->idxCol = (int*)malloc((numNzAIneq + numNzASoc) * sizeof(int));
- socpInfo->ACone->v = (double*)malloc((numNzAIneq + numNzASoc) * sizeof(double));
- socpInfo->ACone->m = numRowAIneq + numRowASoc;
- socpInfo->ACone->n = numOptVar;
- socpInfo->AEqua1 = (cmscp_ccsSpmat*)malloc(sizeof(cmscp_ccsSpmat));
- socpInfo->ACone1 = (cmscp_ccsSpmat*)malloc(sizeof(cmscp_ccsSpmat));
- socpInfo->AEqua1->idxRow = (int*)malloc(numNzAEqua * sizeof(int));
- socpInfo->AEqua1->ptrCol = (int*)malloc((numOptVar+1) * sizeof(int));
- socpInfo->AEqua1->ptrCol[numOptVar] = numNzAEqua;
- socpInfo->AEqua1->v = (double*)malloc(numNzAEqua * sizeof(double));
- socpInfo->AEqua1->m = numRowAEqua;
- socpInfo->AEqua1->n = numOptVar;
- socpInfo->ACone1->idxRow = (int*)malloc((numNzAIneq + numNzASoc) * sizeof(int));
- socpInfo->ACone1->ptrCol = (int*)malloc( (numOptVar + 1) * sizeof(int));
- socpInfo->ACone1->ptrCol[numOptVar] = numNzAIneq + numNzASoc;
- socpInfo->ACone1->v = (double*)malloc((numNzAIneq + numNzASoc) * sizeof(double));
- socpInfo->ACone1->m = numRowAIneq + numRowASoc;
- socpInfo->ACone1->n = numOptVar;
- socpInfo->dimCone = (cmscp_dimCone*)malloc(sizeof(cmscp_dimCone));
- socpInfo->dimCone->l = numRowAIneq;
- socpInfo->dimCone->nsoc = numSubASoc;
- socpInfo->dimCone->q = (int*)malloc(numSubASoc * sizeof(int));
- socpInfo->numOptVar = numOptVar;
- socpInfo->numPhaseOptVar = numPhaseOptVar;
- socpInfo->c = c;
- socpInfo->cntNzAEqua = 0;
- socpInfo->cntRowAEqua = 0;
- socpInfo->cntNzAIneq = 0;
- socpInfo->cntRowAIneq = 0;
- socpInfo->cntNzASoc = numNzAIneq;
- socpInfo->cntRowASoc = numRowAIneq;
- socpInfo->cntSubASoc = 0;
- // 输出
- *socpInfo0 = socpInfo;
- }
- void socpInfoSparseObjEvent(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_sparseInfo* sparseInfo, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup)
- {
- // 函数功能:获取与目标函数和事件约束相关的二阶锥问题参数
- // 输入:
- // socpInfo0:无获目标函数和事件约束相关的二阶锥问题参数
- // refTraj:参考轨迹
- // consTypeInfo:约束类型信息
- // sparseInfo:稀疏(非零元位置)信息
- // setup:问题设置参数
- // 输出:
- // socpInfo0:有目标函数和事件约束相关的二阶锥问题参数
- // 假设函数的自变量按照[x, u, p]排序
- // 假设优化变量的排列顺序为
- // [{x1}, { u1 }, t01, tf1, p1, Tr1, .., { xi }, { ui }, t0i, tfi, pi, Tri, ..., { xN }, { uN }, t0N, tfN, pN, TrN];
- // 其中{ xi } = { xi_1,...,xi_j,...,xi_n }, xi_j表示第i段的第j个状态向量
- // { ui } = { ui_1,...,ui_j,...,ui_n }, xi_j表示第i段的第j个控制输入
- // t0i, tfi表示第i段初始和终端时刻
- // pi表示第i段参数变量
- // Tri表示软信赖域边界
- // 谢磊:2022 / 6 / 27编写
- // 变量声明
- int numPhase, numOptVar, * numPhaseOptVar, iphase, numMesh, numState, numControl, numParameter, numRowAEqua, numNzAEqua, numRowAIneq, numNzAIneq, numRowASoc, numNzASoc, numSubASoc, numNlPath, numSubSocPath, * dimSocPath, * contStateNzPtrRow, * contControlNzPtrRow, * contParameterNzPtrRow, numNzPhaseAEqua, numRowPhaseAIneq, numRowPhaseAEqua, numNzPhaseAIneq, iType, numType, * idxType, cnt, idxRow, numPhaseSubASoc, numRowPointASoc, numRowPhaseASoc, idxStartRowSoc, idxEndRowSoc, numNzPhaseASoc, numNlEndp, numSubSocEndp, * dimSocEndp, * endpIniStateNzPtrRow, * endpFinStateNzPtrRow, * endpIniTimeNzPtrRow, * endpFinTimeNzPtrRow, * endpParameterNzPtrRow, numNlLink, numSubSocLink, * dimSocLink, * linkLeftStatefNzPtrRow, * linkLeftTfNzPtrRow, * linkLeftParameterNzPtrRow, * linkRightState0NzPtrRow, * linkRightT0NzPtrRow, * linkRightParameterNzPtrRow, numNlEvent, numSubSocEvent, * dimSocEvent, idxStartRowSocEvent, * objEventIniStateNzPos, * objEventFinStateNzPos, * objEventIniTimeNzPos, * objEventFinTimeNzPos, * ptrRowAEqua, * idxColAEqua, * ptrRowACone, * idxColACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, cntSubASoc, idxOptVarPhaseStart, cntRowGrd, * iniStateNzPtrRow, * finStateNzPtrRow, * iniTimeNzPtrRow, * finTimeNzPtrRow, * parameterNzPtrRow, * iniStateNzIdxCol, * finStateNzIdxCol, * iniTimeNzIdxCol, * finTimeNzIdxCol, * parameterNzIdxCol, idxOptVarState0, idxOptVarStatef, idxOptVarT0, idxOptVarTf, idxOptVarParameter, idxOptVarTrust, j, nnzType1, i, i0, cntRowAIneqLower, numType1, idxGrdCol, p, numSocEvent, numRowEvent, numRowSocEvent;
- double* c, * valueAEqua, * bEqua, * valueACone, * bCone, w, * eventNlLower, * eventNlUpper, * refState0, * refStatef, refT0, refTf, * refParameter, temp, * eventGrd, * eventRef, *eventSocUpper;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_bound* bounds;
- cmscp_trust* trust;
- cmscp_consType* consType;
- cmscp_nzPos_endp** objEventNzPos;
- cmscp_derInfo* objDerInfo;
- cmscp_derInfo* eventDerInfo;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- // 1. 基本参数
- dimension = setup->dimension;// 问题维数信息
- mesh = setup->mesh;// 网格信息
- numPhase = setup->numPhase;// 阶段数
- bounds = setup->bounds;// 边界信息
- trust = setup->trust;// 信赖域信息
- socpInfo = (*socpInfo0);
- numPhaseOptVar = socpInfo->numPhaseOptVar;// 每段优化变量个数
- objEventNzPos = sparseInfo->objEventNzPos;// 稀疏信息
- consType = consTypeInfo->event;// 约束类型信息
- // 二阶锥规划参数
- numOptVar = socpInfo->numOptVar;
- c = socpInfo->c;
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- cntSubASoc = socpInfo->cntSubASoc;
- //1.1. 梯度信息
- gradInfoSparseObjEvent(&objDerInfo, &eventDerInfo, sparseInfo, refTraj, setup);
- eventGrd = eventDerInfo->der;
- eventRef = eventDerInfo->ref;
- //// 2->目标函数
- floatVecFillin(c, numOptVar, 0);
- // 索引
- idxOptVarPhaseStart = 0; // 各段起始索引
- cntRowGrd = 0; // 目标函数梯度计数
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 各段维数信息
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- if (strcmp(setup->trust->model, "hard") == 0)
- {
- w = 0;
- }
- else if ((strcmp(setup->trust->model, "soft") == 0))
- {
- w = trust->phase[iphase]->softPenaWeight;
- }
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- iniTimeNzIdxCol = objEventNzPos[iphase]->iniTime->idxCol;
- finTimeNzIdxCol = objEventNzPos[iphase]->finTime->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- idxOptVarTrust = idxOptVarParameter + numParameter;
- // x0 前面系数
- for (cnt = iniStateNzPtrRow[0]; cnt < iniStateNzPtrRow[1]; cnt++)
- {
- j = iniStateNzIdxCol[cnt];
- c[idxOptVarState0 + j] = objDerInfo->der[cntRowGrd + j];
- }
- cntRowGrd = cntRowGrd + numState;
- // xf 前面系数
- for (cnt = finStateNzPtrRow[0]; cnt < finStateNzPtrRow[1]; cnt++)
- {
- j = finStateNzIdxCol[cnt];
- c[idxOptVarStatef + j] = objDerInfo->der[cntRowGrd + j];
- }
- cntRowGrd = cntRowGrd + numState;
- // t0 前面系数
- for (cnt = iniTimeNzPtrRow[0]; cnt < iniTimeNzPtrRow[1]; cnt++)
- {
- j = iniTimeNzIdxCol[cnt];
- c[idxOptVarT0] = objDerInfo->der[cntRowGrd + j];
- }
- cntRowGrd = cntRowGrd + 1;
- // tf 前面系数
- for (cnt = finTimeNzPtrRow[0]; cnt < finTimeNzPtrRow[1]; cnt++)
- {
- j = finTimeNzIdxCol[cnt];
- c[idxOptVarTf] = objDerInfo->der[cntRowGrd + j];
- }
- cntRowGrd = cntRowGrd + 1;
- // p 前面系数
- for (cnt = parameterNzPtrRow[0]; cnt < parameterNzPtrRow[1]; cnt++)
- {
- j = parameterNzIdxCol[cnt];
- c[idxOptVarParameter + j] = objDerInfo->der[cntRowGrd + j];
- }
- cntRowGrd = cntRowGrd + numParameter;
- // Tr 前面系数
- c[idxOptVarTrust] = w;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- //// 3-> 事件函数
- // 3->1 类型1:双不等式约束项
- // 上界: A * x0 + B * xf + C * t0 + D * tf + F * p <= b
- // 下界: - A * x0 - B * xf - C * t0 - D * tf - F * p <= b
- iType = 1;
- numType = consType->type[iType].num;
- idxType = consType->type[iType].idx;
- // 维数信息
- numNlEvent = dimension->event->nl;
- numSocEvent = dimension->event->soc->n;
- dimSocEvent = dimension->event->soc->q;
- if (numSocEvent)
- {
- numRowSocEvent = cmscp_sum(dimSocEvent, numSocEvent);
- numRowEvent = numNlEvent + numRowSocEvent;
- }
- else
- {
- numRowEvent = numNlEvent;
- }
- // 约束边界
- eventNlLower = bounds->event->nl->lower;
- eventNlUpper = bounds->event->nl->upper;
- eventSocUpper = bounds->event->soc;
- if (numType)
- {
- nnzType1 = getType1EventNnz(consType->type[iType], objEventNzPos, numPhase); // 第1类约束的一半非零元个数
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt]; // 在Event矩阵的行
- i0 = i + 1; // 在objEvent矩阵的行(Event第一行前加obj)
- cntRowAIneq = cntRowAIneq + 1;
- cntRowAIneqLower = cntRowAIneq + numType; // 下界部分的行位置
- idxOptVarPhaseStart = 0; // 各段起始索引
- idxGrdCol = 0; // 雅可比矩阵的列的计数器
- // 负左端常数项
- temp = -eventRef[i];
- // 每行首元素位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- ptrRowACone[cntRowAIneq+ numType] = cntNzAIneq+ nnzType1;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 离散点个数和各种变量的维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- // 参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // x0 前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p]; // 列坐标
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = eventGrd[i + (idxGrdCol + j) * numRowEvent];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // xf 前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = eventGrd[i + (idxGrdCol + j) * numRowEvent];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // t0 前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarT0;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + 1;
- // tf 前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarTf;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + 1;
- // p 前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = eventGrd[i + (idxGrdCol + j) * numRowEvent];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numParameter;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- // 右端项
- bCone[cntRowAIneq] = temp + eventNlUpper[i]; // 上界
- bCone[numType + cntRowAIneq] = -temp - eventNlLower[i]; // 下界
- cntRowAIneq = cntRowAIneq + 1;
- }
- cntNzAIneq = cntNzAIneq + nnzType1;
- cntRowAIneq = cntRowAIneq + numType;
- }
- // 3->2 类型2:等式约束项
- // A * x0 + B * xf + C * t0 + D * tf + F * p == b
- iType = 1;
- numType = consType->type[iType].num;
- idxType = consType->type[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- i0 = i + 1;
- idxOptVarPhaseStart = 0;
- idxGrdCol = 0; // 雅可比矩阵的列的计数器
- // 负左端常数项
- temp = -eventRef[i];
- // 每行首元素位置
- ptrRowACone[cntRowAEqua] = cntNzAEqua;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 离散点个数和各种变量的维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- // 参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // x0 前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p]; // 列坐标
- cntNzAEqua = cntNzAEqua + 1;
- idxColAEqua[cntNzAEqua] = idxOptVarState0 + j;
- valueACone[cntNzAEqua] = eventGrd[i + (idxGrdCol + j) * numRowEvent];
- temp = temp + valueACone[cntNzAEqua] * refState0[j];
- }
- idxGrdCol = idxGrdCol + numState;
- // xf 前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarStatef + j;
- valueAEqua[cntNzAEqua] = eventGrd[i + (idxGrdCol + j) * numRowEvent];
- temp = temp + valueAEqua[cntNzAEqua] * refStatef[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // t0 前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- idxColAEqua[cntNzAEqua] = idxOptVarT0;
- valueAEqua[cntNzAEqua] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueAEqua[cntNzAEqua] * refT0;
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxGrdCol = idxGrdCol + 1;
- // tf 前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- idxColAEqua[cntNzAEqua] = idxOptVarTf;
- valueAEqua[cntNzAEqua] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueAEqua[cntNzAEqua] * refTf;
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxGrdCol = idxGrdCol + 1;
- // p 前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarParameter + j;
- valueAEqua[cntNzAEqua] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueAEqua[cntNzAEqua] * refParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxGrdCol = idxGrdCol + numParameter;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- // 右端项
- bEqua[cntRowAEqua] = temp + eventNlUpper[i];
- cntRowAEqua = cntRowAEqua + 1;
- }
- }
- // 3->3 类型3:上边界约束项
- // A * x0 + B * xf + C * t0 + D * tf + F * p == b
- iType = 2;
- numType = consType->type[iType].num;
- idxType = consType->type[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- i0 = i + 1;
- idxOptVarPhaseStart = 0;
- idxGrdCol = 0; // 雅可比矩阵的列的计数器
- // 负左端常数项
- temp = -eventRef[i];
- // 每行首元素位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 离散点个数和各种变量的维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- // 参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // x0 前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p]; // 列坐标
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // xf 前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // t0 前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- cntNzAIneq = cntNzAIneq + 1;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- idxColACone[cntNzAIneq] = idxOptVarT0;
- }
- idxGrdCol = idxGrdCol + 1;
- // tf 前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- cntNzAIneq = cntNzAIneq + 1;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- idxColACone[cntNzAIneq] = idxOptVarTf;
- }
- idxGrdCol = idxGrdCol + 1;
- // p 前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numParameter;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- // 右端项
- bEqua[cntRowAIneq] = temp + eventNlUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- // 3->4 类型4:下边界约束项
- // A * x0 + B * xf + C * t0 + D * tf + F * p == b
- iType = 3;
- numType = consType->type[iType].num;
- idxType = consType->type[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- i0 = i + 1;
- idxOptVarPhaseStart = 0;
- idxGrdCol = 0; // 雅可比矩阵的列的计数器
- // 负左端常数项
- temp = -eventRef[i];
- // 每行首元素位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 离散点个数和各种变量的维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- // 参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // x0 前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p]; // 列坐标
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = -eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // xf 前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = -eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numState;
- // t0 前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarT0;
- valueACone[cntNzAIneq] = -eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- }
- idxGrdCol = idxGrdCol + 1;
- // tf 前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarTf;
- valueACone[cntNzAIneq] = -eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + 1;
- // p 前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = -eventGrd[i + idxGrdCol * numRowEvent];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxGrdCol = idxGrdCol + numParameter;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- // 右端项
- bEqua[cntRowAIneq] = temp + eventNlUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- // 3->5 二阶锥约束
- // 事件约束维度信息
- if (numSocEvent)
- {
- for (i = 0; i < numRowSocEvent; i++)
- {
- i0 = numNlEvent + 1 + i;
- temp = -eventRef[numNlEvent + i];
- // 每行首个元素的位置
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxOptVarPhaseStart = 0;// 各段起始索引
- idxGrdCol = 0;// 雅可比矩阵的列的计数
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 离散点个数和各种变量的维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- // 稀疏信息
- iniStateNzPtrRow = objEventNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = objEventNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = objEventNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = objEventNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = objEventNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = objEventNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = objEventNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = objEventNzPos[iphase]->parameter->idxCol;
- // 每段各种所需部分的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarPhaseStart + numMesh * (numState + numControl) + 2;
- // 参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // x0 前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p]; // 列
- cntNzASoc = cntNzASoc + 1;
- idxColACone[cntNzASoc] = idxOptVarState0 + j;
- valueACone[cntNzASoc] = eventGrd[numNlEvent + i + (idxGrdCol + j) * numRowEvent];
- temp = temp + valueACone[cntNzASoc] * refState0[j];
- }
- idxGrdCol = idxGrdCol + numState;
- // xf 前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p]; // 列
- cntNzASoc = cntNzASoc + 1;
- idxColACone[cntNzASoc] = idxOptVarStatef + j;
- valueACone[cntNzASoc] = eventGrd[numNlEvent + i + (idxGrdCol + j) * numRowEvent];
- temp = temp + valueACone[cntNzASoc] * refStatef[j];
- }
- idxGrdCol = idxGrdCol + numState;
- // t0 前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- cntNzASoc = cntNzASoc + 1;
- idxColACone[cntNzASoc] = idxOptVarT0 + 1;
- valueACone[cntNzASoc] = eventGrd[numNlEvent + i + (idxGrdCol) * numRowEvent];
- temp = temp + valueACone[cntNzASoc] * refT0;
- }
- idxGrdCol = idxGrdCol + 1;
- // tf 前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- cntNzASoc = cntNzASoc + 1;
- idxColACone[cntNzASoc] = idxOptVarTf + 1;
- valueACone[cntNzASoc] = eventGrd[numNlEvent + i + (idxGrdCol) * numRowEvent];
- temp = temp + valueACone[cntNzASoc] * refTf;
- }
- idxGrdCol = idxGrdCol + 1;
- // p 前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p]; // 列
- cntNzASoc = cntNzASoc + 1;
- idxColACone[cntNzASoc] = idxOptVarParameter + j;
- valueACone[cntNzASoc] = eventGrd[numNlEvent + i + (idxGrdCol + j) * numRowEvent];
- temp = temp + valueACone[cntNzASoc] * refParameter[j];
- }
- idxGrdCol = idxGrdCol + numParameter;
- // 更新索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- cntRowASoc = cntRowASoc + 1;
- }
- for (i = 0; i < numSocEvent; i++)
- {
- dimCone->q[cntSubASoc] = dimSocEvent[i];
- cntSubASoc = cntSubASoc + 1;
- }
- }
- if (objDerInfo->ref) { free(objDerInfo->ref); }
- if (objDerInfo->der) { free(objDerInfo->der); }
- if (objDerInfo) { free(objDerInfo); }
- if (eventDerInfo->ref) { free(eventDerInfo->ref); }
- if (eventDerInfo->der) { free(eventDerInfo->der); }
- if (eventDerInfo) { free(eventDerInfo); }
- }
- void socpInfoSparseEndp(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_sparseInfo* sparseInfo, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup) {
- // 函数功能:获取端点函数约束的二阶锥规划问题参数
- // 输入:
- // socpInfo0:无获端点函数约束相关的二阶锥问题参数
- // refTraj:参考轨迹
- // consTypeInfo:约束类型信息
- // sparseInfo:稀疏(非零元位置)信息
- // setup:问题设置参数
- // 输出:
- // socpInfo0:有端点函数约束相关的二阶锥问题参数
- // 假设函数的自变量按照[x, u, p]排序
- // 谢磊:2022 / 6 / 27编写
- // 变量声明
- int numPhase, * numPhaseOptVar, iphase, numMesh, numState, numControl, iType, numType, * idxType, cnt, numNlEndp, * dimSocEndp, * idxColAEqua, * ptrRowAEqua, * idxColACone, * ptrRowACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, cntSubASoc, idxOptVarPhaseStart, numSocEndp, numRowSocEndp, * iniStateNzPtrRow, * finStateNzPtrRow, * iniTimeNzPtrRow, * finTimeNzPtrRow, * parameterNzPtrRow, * iniStateNzIdxCol, * finStateNzIdxCol, * parameterNzIdxCol, idxEndpGrdColStatef, idxEndpGrdColT0, idxEndpGrdColTf, idxEndpGrdColParameter, idxOptVarState0, idxOptVarStatef, idxOptVarT0, idxOptVarTf, idxOptVarParameter, nnzType1, i, p, j, numEndpOutput, i0, numOptVar;
- double* valueAEqua, * bEqua, * valueACone, * bCone, * endpNlLower, * endpNlUpper, * endpSocUpper, * refState0, * refStatef, refT0, refTf, * refParameter, * endpRef, * endpGrd, temp, *c;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_bound* bounds;
- cmscp_trust* trust;
- cmscp_consType** consType;
- cmscp_type* consTypePhase;
- cmscp_nzPos_endp** endpNzPos;
- cmscp_derInfo** gradInfo;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- // 1. 基本参数
- dimension = setup->dimension;// 问题维数信息
- mesh = setup->mesh;// 网格信息
- numPhase = setup->numPhase;// 阶段数
- bounds = setup->bounds;// 边界信息
- trust = setup->trust;// 信赖域信息
- socpInfo = (*socpInfo0);
- numOptVar = socpInfo->numOptVar;
- numPhaseOptVar = socpInfo->numPhaseOptVar;// 每段优化变量个数
- endpNzPos = sparseInfo->endpNzPos;// 稀疏信息
- consType = consTypeInfo->endp;// 约束类型信息
- // 二阶锥规划参数
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- cntSubASoc = socpInfo->cntSubASoc;
- // 梯度信息
- gradInfo = (cmscp_derInfo**)malloc(numPhase * sizeof(cmscp_derInfo));
- gradInfoSparseEndp(gradInfo, sparseInfo, refTraj, setup);
-
- //// 填充非零元素
- // 优化变量中各段起始索引
- idxOptVarPhaseStart = 0;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 各段维数信息
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numNlEndp = dimension->phase[iphase]->endpoint->nl;
- numSocEndp = dimension->phase[iphase]->endpoint->soc->n;
- dimSocEndp = dimension->phase[iphase]->endpoint->soc->q;
- if (numSocEndp)
- {
- numRowSocEndp = cmscp_sum(dimSocEndp, numSocEndp);
- }
- else
- {
- numRowSocEndp = 0;
- }
- if (dimension->phase[iphase]->endpoint->soc->n)
- {
- numEndpOutput = numNlEndp + cmscp_sum(dimSocEndp, numSocEndp);
- }
- else
- {
- numEndpOutput = numNlEndp;
- }
- // 各段约束类型信息
- consTypePhase = consType[iphase]->type;
- // 各段稀疏信息
- iniStateNzPtrRow = endpNzPos[iphase]->iniState->ptrRow;
- finStateNzPtrRow = endpNzPos[iphase]->finState->ptrRow;
- iniTimeNzPtrRow = endpNzPos[iphase]->iniTime->ptrRow;
- finTimeNzPtrRow = endpNzPos[iphase]->finTime->ptrRow;
- parameterNzPtrRow = endpNzPos[iphase]->parameter->ptrRow;
- iniStateNzIdxCol = endpNzPos[iphase]->iniState->idxCol;
- finStateNzIdxCol = endpNzPos[iphase]->finState->idxCol;
- parameterNzIdxCol = endpNzPos[iphase]->parameter->idxCol;
- // 各段约束上下界
- endpNlLower = bounds->phase[iphase]->endpoint->nl->lower;
- endpNlUpper = bounds->phase[iphase]->endpoint->nl->upper;
- endpSocUpper = bounds->phase[iphase]->endpoint->soc;
- // 各段参考轨迹
- refState0 = refTraj[iphase]->state->v;
- refStatef = refTraj[iphase]->state->v + (numMesh - 1) * numState;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- // 各段梯度信息
- endpRef = gradInfo[iphase]->ref;
- endpGrd = gradInfo[iphase]->der;
- // 各段所需变量的起始索引
- idxEndpGrdColStatef = numState;// EndpGrd矩阵中终端点状态的对应的列的起始索引
- idxEndpGrdColT0 = 2 * numState;// EndpGrd矩阵中初始时刻的对应的列的起始索引
- idxEndpGrdColTf = idxEndpGrdColT0 + 1;// EndpGrd矩阵中终端时刻的对应的列的起始索引
- idxEndpGrdColParameter = idxEndpGrdColTf + 1;// EndpGrd矩阵中终端时刻的对应的列的起始索引
- idxOptVarState0 = idxOptVarPhaseStart;// 优化变量中每段的初始点状态的索引
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;// 优化变量中每段的终端点状态的索引
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);// 优化变量中每段的初始点时刻的索引
- idxOptVarTf = idxOptVarT0 + 1;// 优化变量中每段的终端点时刻的索引
- idxOptVarParameter = idxOptVarTf + 1;// 优化变量中每段的静态参数的索引
- //// 1 类型1:双不等式约束项
- // 上界:A * x0 + B * xf + C * t0 + D * tf + F * p <= b
- // 下界: - A * x0 - B * xf - C * t0 - D * tf - F * p <= b
- iType = 0;
- numType = consTypePhase[iType].num;
- idxType = consTypePhase[iType].idx;
- if (numType)
- {
- nnzType1 = getType1EndpNnz(consTypePhase[iType], endpNzPos[iphase]);// 一半非零元个数
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- temp = -endpRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- ptrRowACone[cntRowAIneq + numType] = cntNzAIneq + nnzType1;
- // x0前面系数
- for (p = iniStateNzPtrRow[i]; p < iniStateNzPtrRow[i + 1]; p++)
- {
- j = iniStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = endpGrd[i+j * numEndpOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // xf前面系数
- for (p = finStateNzPtrRow[i]; p < finStateNzPtrRow[i + 1]; p++)
- {
- j = finStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = endpGrd[i+(idxEndpGrdColStatef + j) * numEndpOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // t0前面系数
- if (iniTimeNzPtrRow[i] + 1 == iniTimeNzPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarT0;
- valueACone[cntNzAIneq] = endpGrd[i+idxEndpGrdColT0 * numEndpOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // tf前面系数
- if (finTimeNzPtrRow[i] + 1 == finTimeNzPtrRow[i + 1]) {
- idxColACone[cntNzAIneq] = idxOptVarTf;
- valueACone[cntNzAIneq] = endpGrd[i+idxEndpGrdColTf * numEndpOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i]; p < parameterNzPtrRow[i + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = endpGrd[i+(idxEndpGrdColParameter + j) * numEndpOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp + endpNlUpper[i];
- bCone[cntRowAIneq + numType] = -temp - endpNlLower[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- cntNzAIneq = cntNzAIneq + nnzType1;
- cntRowAIneq = cntRowAIneq + numType;
- }
- //// 2 类型2:等式约束项,表达式:
- // A * x0 + B * xf + C * t0 + D * tf + F * p == b
- iType = 1;
- numType = consTypePhase[iType].num;
- idxType = consTypePhase[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- temp = -endpRef[i];
- // 每行首个非零元的位置
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- // x0前面系数
- for (p = iniStateNzPtrRow[i]; p < iniStateNzPtrRow[i + 1]; p++)
- {
- j = iniStateNzIdxCol[p];// 列
- idxColAEqua[cntNzAEqua] = idxOptVarState0 + j;
- valueAEqua[cntNzAEqua] = endpGrd[i+j * numEndpOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refState0[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // xf前面系数
- for (p = finStateNzPtrRow[i]; p < finStateNzPtrRow[i + 1]; p++)
- {
- j = finStateNzIdxCol[p];// 列
- idxColAEqua[cntNzAEqua] = idxOptVarStatef + j;
- valueAEqua[cntNzAEqua] = endpGrd[i+(idxEndpGrdColStatef + j) * numEndpOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refStatef[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // t0前面系数
- if (iniTimeNzPtrRow[i] + 1 == iniTimeNzPtrRow[i + 1])
- {
- idxColAEqua[cntNzAEqua] = idxOptVarT0;
- valueAEqua[cntNzAEqua] = endpGrd[i+idxEndpGrdColT0 * numEndpOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refT0;
- cntNzAEqua = cntNzAEqua + 1;
- }
- // tf前面系数
- if (finTimeNzPtrRow[i] + 1 == finTimeNzPtrRow[i + 1])
- {
- idxColAEqua[cntNzAEqua] = idxOptVarTf;
- valueAEqua[cntNzAEqua] = endpGrd[i+idxEndpGrdColTf * numEndpOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refTf;
- cntNzAEqua = cntNzAEqua + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i]; p < parameterNzPtrRow[i + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxColAEqua[cntNzAEqua] = idxOptVarParameter + j;
- valueAEqua[cntNzAEqua] = endpGrd[i+(idxEndpGrdColParameter + j) * numEndpOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- bEqua[cntRowAEqua] = temp + endpNlUpper[i];
- cntRowAEqua = cntRowAEqua + 1;
- }
- }
- //// 类型3:上不等式约束项,表达式:
- // A* x0 + B * xf + C * t0 + D * tf + F * p <= b
- iType = 2;
- numType = consTypePhase[iType].num;
- idxType = consTypePhase[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- temp = -endpRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // x0前面系数
- for (p = iniStateNzPtrRow[i]; p < iniStateNzPtrRow[i + 1]; p++)
- {
- j = iniStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = endpGrd[i+j * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // xf前面系数
- for (p = finStateNzPtrRow[i]; p < finStateNzPtrRow[i + 1]; p++)
- {
- j = finStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = endpGrd[i+(idxEndpGrdColStatef + j) * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // t0前面系数
- if (iniTimeNzPtrRow[i] + 1 == iniTimeNzPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarT0;
- valueACone[cntNzAIneq] = endpGrd[i+idxEndpGrdColT0 * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // tf前面系数
- if (finTimeNzPtrRow[i] + 1 == finTimeNzPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarTf;
- valueACone[cntNzAIneq] = endpGrd[i+idxEndpGrdColTf * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i]; p < parameterNzPtrRow[i + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = endpGrd[i+(idxEndpGrdColParameter + j) * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp + endpNlUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- //// 类型4:下不等式约束项
- // 表达式: - A * x0 - B * xf - C * t0 - D * tf - F * p <= -b
- iType = 3;
- numType = consTypePhase[iType].num;
- idxType = consTypePhase[iType].idx;
- if (numType)
- {
- for (cnt = 0; cnt < numType; cnt++)
- {
- i = idxType[cnt];
- temp = endpRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // x0前面系数
- for (p = iniStateNzPtrRow[i]; p < iniStateNzPtrRow[i + 1]; p++)
- {
- j = iniStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarState0 + j;
- valueACone[cntNzAIneq] = -endpGrd[i+j * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // xf前面系数
- for (p = finStateNzPtrRow[i]; p < finStateNzPtrRow[i + 1]; p++)
- {
- j = finStateNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarStatef + j;
- valueACone[cntNzAIneq] = -endpGrd[i+(idxEndpGrdColStatef + j) * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // t0前面系数
- if (iniTimeNzPtrRow[i] + 1 == iniTimeNzPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarT0;
- valueACone[cntNzAIneq] = -endpGrd[i+idxEndpGrdColT0 * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // tf前面系数
- if (finTimeNzPtrRow[i] + 1 == finTimeNzPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarTf;
- valueACone[cntNzAIneq] = -endpGrd[i+idxEndpGrdColTf * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i]; p < parameterNzPtrRow[i + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = -endpGrd[i+(idxEndpGrdColParameter + j) * numEndpOutput];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp - endpNlLower[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- if (numSocEndp)
- {
- for (i = 0; i < numRowSocEndp; i++)
- {
- i0 = numNlEndp + i;
- temp = -endpRef[i0];
- // 每行首个非零元的位置
- ptrRowACone[cntRowASoc] = cntNzASoc;
- // x0前面系数
- for (p = iniStateNzPtrRow[i0]; p < iniStateNzPtrRow[i0 + 1]; p++)
- {
- j = iniStateNzIdxCol[p];// 列
- idxColACone[cntNzASoc] = idxOptVarState0 + j;
- valueACone[cntNzASoc] = endpGrd[i0+j * numEndpOutput];
- temp = temp + valueACone[cntNzASoc] * refState0[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // xf前面系数
- for (p = finStateNzPtrRow[i0]; p < finStateNzPtrRow[i0 + 1]; p++)
- {
- j = finStateNzIdxCol[p];// 列
- idxColACone[cntNzASoc] = idxOptVarStatef + j;
- valueACone[cntNzASoc] = endpGrd[i0+(idxEndpGrdColStatef + j) * numEndpOutput];
- temp = temp + valueACone[cntNzASoc] * refStatef[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // t0前面系数
- if (iniTimeNzPtrRow[i0] + 1 == iniTimeNzPtrRow[i0 + 1])
- {
- idxColACone[cntNzASoc] = idxOptVarT0;
- valueACone[cntNzASoc] = endpGrd[i0+idxEndpGrdColT0 * numEndpOutput];
- temp = temp + valueACone[cntNzASoc] * refT0;
- cntNzASoc = cntNzASoc + 1;
- }
- // tf前面系数
- if (finTimeNzPtrRow[i0] + 1 == finTimeNzPtrRow[i0 + 1])
- {
- idxColACone[cntNzASoc] = idxOptVarTf;
- valueACone[cntNzASoc] = endpGrd[i0+idxEndpGrdColTf * numEndpOutput];
- temp = temp + valueACone[cntNzASoc] * refTf;
- cntNzASoc = cntNzASoc + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxColACone[cntNzASoc] = idxOptVarParameter + j;
- valueACone[cntNzASoc] = endpGrd[i0+(idxEndpGrdColParameter + j) * numEndpOutput];
- temp = temp + valueACone[cntNzASoc] * refParameter[j];
- cntNzASoc = cntNzASoc + 1;
- }
- bCone[cntRowASoc] = temp + endpSocUpper[i];
- cntRowASoc = cntRowASoc + 1;
- }
- for (i = 0; i < numSocEndp; i++)
- {
- dimCone->q[cntSubASoc] = dimSocEndp[i];
- cntSubASoc = cntSubASoc + 1;
- }
- }
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- free(gradInfo[iphase]->ref);
- free(gradInfo[iphase]->der);
- free(gradInfo[iphase]);
- }
- socpInfo->cntNzAEqua = cntNzAEqua;
- socpInfo->cntRowAEqua = cntRowAEqua;
- socpInfo->cntNzAIneq = cntNzAIneq;
- socpInfo->cntRowAIneq = cntRowAIneq;
- socpInfo->cntNzASoc = cntNzASoc;
- socpInfo->cntRowASoc = cntRowASoc;
- socpInfo->cntSubASoc = cntSubASoc;
- free(gradInfo);
- }
- void socpInfoSparseCont(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_sparseInfo* sparseInfo, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup)
- {
- // 函数功能:获取连续函数约束的二阶锥规划问题参数
- // 输入:
- // socpInfo:二阶锥信息结构体
- // refTraj:参考轨迹信息
- // consTypeInfo:约束类型信息
- // sparseInfo:稀疏(非零元位置)信息
- // setup:问题设置参数
- // 输出:
- // socpInfo:二阶锥信息结构体
- // 谢磊:2022 / 4 / 18编写
- // 变量声明
- int numPhase, * numPhaseOptVar, iphase, numMesh, numState, numControl, cnt, * idxColAEqua, * ptrRowAEqua, * idxColACone, * ptrRowACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, cntSubASoc, idxOptVarPhaseStart, * parameterNzPtrRow, * parameterNzIdxCol, idxOptVarT0, idxOptVarTf, idxOptVarParameter, nnzType1, i, p, j, i0, numParameter, numNlPath, numSocPath, * dimSocPath, numRowSocPath, numRowPath, * stateNzPtrRow, * controlNzPtrRow, * stateNzIdxCol, * controlNzIdxCol, * stateNzIdxRow, * controlNzIdxRow, * parameterNzIdxRow, * stateNzPtrCol, * controlNzPtrCol, * parameterNzPtrCol, numType1, numType2, numType3, numType4, * idxType1, * idxType2, * idxType3, * idxType4, idxPathGrdColControl, idxPathGrdColParameter, idxDynGrdColControl, idxDynGrdColParameter, idxOptVarStateK, idxOptVarControlK, k, idxDynGrd, idxOptVarStateKm1, idxOptVarControlKm1, numContInput, idxPathGrd;;
- double* valueAEqua, * bEqua, * valueACone, * bCone, refT0, refTf, * refParameter, temp, * refState, * refControl, * meshNode, * refStateK, * refControlK, * pathNlLower, * pathNlUpper, * pathSocUpper, * dynRefK, * dynGrdK, * pathRefK, * pathGrdK, * dynRefKm1, * dynGrdKm1, * pathRefKm1, * pathGrdKm1, a1, a2, * refStateKm1, * refControlKm1, * tempPtr;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_bound* bounds;
- cmscp_consType** consType;
- cmscp_nzPos_cont** contNzPos;
- cmscp_derInfo** dynGradInfo;
- cmscp_derInfo** pathGradInfo;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- // 分配空间
- dynGradInfo = (cmscp_derInfo**)malloc(sizeof(cmscp_derInfo*));
- dynGradInfo[0] = (cmscp_derInfo*)malloc(sizeof(cmscp_derInfo));
- pathGradInfo = (cmscp_derInfo**)malloc(sizeof(cmscp_derInfo*));
- pathGradInfo[0] = (cmscp_derInfo*)malloc(sizeof(cmscp_derInfo));
- // //基本参数
- dimension = setup->dimension;// 问题维数信息
- mesh = setup->mesh;// 网格信息
- bounds = setup->bounds;// 边界信息
- numPhase = setup->numPhase;// 阶段数
- socpInfo = *socpInfo0;
- numPhaseOptVar = socpInfo->numPhaseOptVar;// 每段优化变量个数
- contNzPos = sparseInfo->contNzPos;// 稀疏信息
- consType = consTypeInfo->path;// 约束类型信息
- // 二阶锥规划参数
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- cntSubASoc = socpInfo->cntSubASoc;
- // 优化变量中各段起始索引
- idxOptVarPhaseStart = 0;
- for (iphase = 0; iphase < numPhase; iphase++)
- {
- // 维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- numNlPath = dimension->phase[iphase]->path->nl;
- numSocPath = dimension->phase[iphase]->path->soc->n;
- dimSocPath = dimension->phase[iphase]->path->soc->q;
- if (numSocPath)
- {
- numRowSocPath = cmscp_sum(dimSocPath, numSocPath);
- }
- else
- {
- numRowSocPath = 0;
- }
- numRowPath = numNlPath + numRowSocPath;
- // 约束上下界简写
- pathNlLower = bounds->phase[iphase]->path->nl->lower;
- pathNlUpper = bounds->phase[iphase]->path->nl->upper;
- pathSocUpper = bounds->phase[iphase]->path->soc;
- // 参考轨迹
- refState = refTraj[iphase]->state->v;
- refControl = refTraj[iphase]->control->v;
- refParameter = refTraj[iphase]->parameter;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- // 稀疏信息
- stateNzPtrRow = contNzPos[iphase]->state->ptrRow;
- controlNzPtrRow = contNzPos[iphase]->control->ptrRow;
- parameterNzPtrRow = contNzPos[iphase]->parameter->ptrRow;
- stateNzIdxCol = contNzPos[iphase]->state->idxCol;
- controlNzIdxCol = contNzPos[iphase]->control->idxCol;
- parameterNzIdxCol = contNzPos[iphase]->parameter->idxCol;
- stateNzIdxRow = contNzPos[iphase]->state->idxRow;
- controlNzIdxRow = contNzPos[iphase]->control->idxRow;
- parameterNzIdxRow = contNzPos[iphase]->parameter->idxRow;
- stateNzPtrCol = contNzPos[iphase]->state->ptrCol;
- controlNzPtrCol = contNzPos[iphase]->control->ptrCol;
- parameterNzPtrCol = contNzPos[iphase]->parameter->ptrCol;
- // 约束类型简写
- numType1 = consType[iphase]->type[0].num;
- numType2 = consType[iphase]->type[1].num;
- numType3 = consType[iphase]->type[2].num;
- numType4 = consType[iphase]->type[3].num;
- idxType1 = consType[iphase]->type[0].idx;
- idxType2 = consType[iphase]->type[1].idx;
- idxType3 = consType[iphase]->type[2].idx;
- idxType4 = consType[iphase]->type[3].idx;
- nnzType1 = getType1PathNnz(consType[iphase]->type[0], contNzPos[iphase], numState); // Type1约束的一半非零元个数
- //离散节点
- meshNode = mesh[iphase]->node;
- // 索引
- idxPathGrdColControl = numRowPath * numState;
- idxPathGrdColParameter = numRowPath * (numState + numControl);
- idxDynGrdColControl = numState * numState;
- idxDynGrdColParameter = numState * (numState + numControl);
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarTf + 1;
- idxOptVarStateK = idxOptVarPhaseStart;
- idxOptVarControlK = idxOptVarPhaseStart + numMesh * numState;
- // 导数信息空间
- numContInput = numState + numControl + numParameter;
- dynRefKm1 = (double*)malloc(numState * sizeof(double));
- dynGrdKm1 = (double*)malloc(numState * numContInput * sizeof(double));
- pathRefKm1 = (double*)malloc(numRowPath * sizeof(double));
- pathGrdKm1 = (double*)malloc(numRowPath * numContInput * sizeof(double));
- dynRefK = (double*)malloc(numState * sizeof(double));
- dynGrdK = (double*)malloc(numState * numContInput * sizeof(double));
- pathRefK = (double*)malloc(numRowPath * sizeof(double));
- pathGrdK = (double*)malloc(numRowPath * numContInput * sizeof(double));
- for (k = 0; k < numMesh; k++)
- {
- refStateK = refState + k * numState;
- refControlK = refControl + k * numControl;
- refStateKm1 = refState + (k - 1) * numState;
- refControlKm1 = refControl + (k - 1) * numControl;
- // 梯度信息
- dynGradInfo[0]->ref = dynRefK;
- dynGradInfo[0]->der = dynGrdK;
- pathGradInfo[0]->ref = pathRefK;
- pathGradInfo[0]->der = pathGrdK;
- gradInfoSparseCont(dynGradInfo, pathGradInfo, refStateK, refControlK,
- refParameter, stateNzIdxRow, controlNzIdxRow,
- parameterNzIdxRow, stateNzPtrCol, controlNzPtrCol,
- parameterNzPtrCol, numState,
- numControl, numParameter, numRowPath,
- iphase, setup);
- if (k > 0) {
- // 常数项
- a1 = 0.5 * (meshNode[k] - meshNode[k - 1]);
- a2 = a1 * (refTf - refT0);
- for (i = 0; i < numState; i++) { // 行
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- temp = 0;
- // xK前面系数
- for (p = stateNzPtrRow[i]; p < stateNzPtrRow[i + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxDynGrd = j * numState + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarStateKm1 + j;
- valueAEqua[cntNzAEqua] = a2 * dynGrdKm1[idxDynGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refStateKm1[j];
- if (i == j)
- {
- valueAEqua[cntNzAEqua] = valueAEqua[cntNzAEqua] + 1;
- }
- cntNzAEqua = cntNzAEqua + 1;
- }
- // xKp1前面系数
- for (p = stateNzPtrRow[i]; p < stateNzPtrRow[i + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxDynGrd = j * numState + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarStateK + j;
- valueAEqua[cntNzAEqua] = a2 * dynGrdK[idxDynGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refStateK[j];
- if (i == j)
- {
- valueAEqua[cntNzAEqua] = valueAEqua[cntNzAEqua] - 1;
- }
- cntNzAEqua = cntNzAEqua + 1;
- }
- // uK前面系数
- for (p = controlNzPtrRow[i]; p < controlNzPtrRow[i + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxDynGrd = idxDynGrdColControl + j * numState + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarControlKm1 + j;
- valueAEqua[cntNzAEqua] = a2 * dynGrdKm1[idxDynGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refControlKm1[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // uKp1前面系数
- for (p = controlNzPtrRow[i]; p < controlNzPtrRow[i + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxDynGrd = idxDynGrdColControl + j * numState + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarControlK + j;
- valueAEqua[cntNzAEqua] = a2 * dynGrdK[idxDynGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refControlK[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // t0前面系数
- idxColAEqua[cntNzAEqua] = idxOptVarT0;
- valueAEqua[cntNzAEqua] = -a1 * (dynRefKm1[i] + dynRefK[i]);
- cntNzAEqua = cntNzAEqua + 1;
- // tf前面系数
- idxColAEqua[cntNzAEqua] = idxOptVarTf;
- valueAEqua[cntNzAEqua] = -valueAEqua[cntNzAEqua - 1];
- cntNzAEqua = cntNzAEqua + 1;
- // p前面系数
- for (p = parameterNzPtrRow[i]; p < parameterNzPtrRow[i + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxDynGrd = idxDynGrdColParameter + j * numState + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarParameter + j;
- valueAEqua[cntNzAEqua] = a2 * (dynGrdKm1[idxDynGrd] + dynGrdK[idxDynGrd]);
- temp = temp + valueAEqua[cntNzAEqua] * refParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- bEqua[cntRowAEqua] = temp;
- cntRowAEqua = cntRowAEqua + 1;
- }
- }
- //// 路径约束的系数和右端项
- // 类型1:双不等式约束项
- // 上界: Ak* xK + BK * uK + CK * p <= bK
- // 下界: - Ak * xK - BK * uK - CK * p <= -bK
-
- if (numType1)
- {
- for (cnt = 0; cnt < numType1; cnt++)
- {
- i = idxType1[cnt];// 在路径约束函数中的行位置
- i0 = i + numState;// 在连续约束函数中的行位置
- temp = -pathRefK[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- ptrRowACone[cntRowAIneq + numType1] = cntNzAIneq + nnzType1;
- // xK前面系数矩阵
- for (p = stateNzPtrRow[i0]; p < stateNzPtrRow[i0 + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxPathGrd = j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarStateK + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refStateK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // uK前面系数矩阵
- for (p = controlNzPtrRow[i0]; p < controlNzPtrRow[i0 + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColControl + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarControlK + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refControlK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数矩阵
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++) {
- j = parameterNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColParameter + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右端项
- bCone[cntRowAIneq] = temp + pathNlUpper[i]; // 上界
- bCone[numType1 + cntRowAIneq] = -temp - pathNlLower[i]; // 下界
- cntRowAIneq = cntRowAIneq + 1;
- }
- cntNzAIneq = cntNzAIneq + nnzType1;
- cntRowAIneq = cntRowAIneq + numType1;
- }
- // 等式约束项,表达式:
- // Ak* xK + BK * uK + CK * p == bK
- if (numType2)
- {
- for (cnt = 0; cnt < numType2; cnt++)
- {
- i = idxType2[cnt];
- i0 = i + numState;
- temp = -pathRefK[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- // xK前面系数
- for (p = stateNzPtrRow[i0]; p < stateNzPtrRow[i0 + 1]; p++) {
- j = stateNzIdxCol[p];// 列
- idxPathGrd = j * numRowPath + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarStateK + j;
- valueAEqua[cntNzAEqua] = pathGrdK[idxPathGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refStateK[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // uK前面系数
- for (p = controlNzPtrRow[i0]; p < controlNzPtrRow[i0 + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColControl + j * numRowPath + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarControlK + j;
- valueAEqua[cntNzAEqua] = pathGrdK[idxPathGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refControlK[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColParameter + j * numRowPath + i;// 第(i, j)个元素
- idxColAEqua[cntNzAEqua] = idxOptVarParameter + j;
- valueAEqua[cntNzAEqua] = pathGrdK[idxPathGrd];
- temp = temp + valueAEqua[cntNzAEqua] * refParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- bEqua[cntRowAEqua] = temp + pathNlUpper[i];
- cntRowAEqua = cntRowAEqua + 1;
- }
- }
- // 上单边不等式约束项,表达式:
- // Ak* xK + BK * uK + CK * p <= bK
- if (numType3)
- {
- for (cnt = 0; cnt < numType3; cnt++)
- {
- i = idxType3[cnt];
- i0 = i + numState;
- temp = -pathRefK[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // xK前面系数
- for (p = stateNzPtrRow[i0]; p < stateNzPtrRow[i0 + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxPathGrd = j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarStateK + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refStateK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // uK前面系数
- for (p = controlNzPtrRow[i0]; p < controlNzPtrRow[i0 + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColControl + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarControlK + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refControlK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColParameter + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp + pathNlUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- // 下单边不等式约束项,表达式:
- // -Ak * xK - BK * uK - CK * p <= -bK
- if (numType4) {
- for (cnt = 0; cnt < numType4; cnt++)
- {
- i = idxType4[cnt];
- i0 = i + numState;
- temp = pathRefK[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // xK前面系数
- for (p = stateNzPtrRow[i0]; p < stateNzPtrRow[i0 + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxPathGrd = j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarStateK + j;
- valueACone[cntNzAIneq] = -pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refStateK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // uK前面系数
- for (p = controlNzPtrRow[i0]; p < controlNzPtrRow[i0 + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColControl + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarControlK + j;
- valueACone[cntNzAIneq] = -pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refControlK[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColParameter + j * numRowPath + i;// 第(i, j)个元素
- idxColACone[cntNzAIneq] = idxOptVarParameter + j;
- valueACone[cntNzAIneq] = -pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzAIneq] * refParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp - pathNlLower[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- if (numSocPath)
- {
- for (i = 0; i < numRowSocPath; i++)
- {
- i0 = numState + numNlPath + i;
- temp = -pathRefK[numNlPath + i];
- ptrRowACone[cntRowASoc] = cntNzASoc;
- // xK前面系数
- for (p = stateNzPtrRow[i0]; p < stateNzPtrRow[i0 + 1]; p++)
- {
- j = stateNzIdxCol[p];// 列
- idxPathGrd = j * numRowPath + numNlPath + i;// 第(i, j)个元素
- idxColACone[cntNzASoc] = idxOptVarStateK + j;
- valueACone[cntNzASoc] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzASoc] * refStateK[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // uK前面系数
- for (p = controlNzPtrRow[i0]; p < controlNzPtrRow[i0 + 1]; p++)
- {
- j = controlNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColControl + j * numRowPath + numNlPath + i;// 第(i, j)个元素
- idxColACone[cntNzASoc] = idxOptVarControlK + j;
- valueACone[cntNzASoc] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzASoc] * refControlK[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // p前面系数
- for (p = parameterNzPtrRow[i0]; p < parameterNzPtrRow[i0 + 1]; p++)
- {
- j = parameterNzIdxCol[p];// 列
- idxPathGrd = idxPathGrdColParameter + j * numRowPath + numNlPath + i;// 第(i, j)个元素
- idxColACone[cntNzASoc] = idxOptVarParameter + j;
- valueACone[cntNzASoc] = pathGrdK[idxPathGrd];
- temp = temp + valueACone[cntNzASoc] * refParameter[j];
- cntNzASoc = cntNzASoc + 1;
- }
- bCone[cntRowASoc] = temp + pathSocUpper[i];
- cntRowASoc = cntRowASoc + 1;
- }
- for (i = 0; i < numSocPath; i++)
- {
- cntSubASoc = cntSubASoc + 1;
- dimCone->q[cntSubASoc] = dimSocPath[i];
- }
- }
- // dynRefKm1和dynRefK空间互换,实现的功能为dynRefKm1=dynRefK,且不需要另外创建空间;
- tempPtr = dynRefKm1; dynRefKm1 = dynRefK; dynRefK = tempPtr;
- tempPtr = dynGrdKm1; dynGrdKm1 = dynGrdK; dynGrdK = tempPtr;
- tempPtr = pathRefKm1; pathRefKm1 = pathRefK; pathRefK = tempPtr;
- tempPtr = pathGrdKm1; pathGrdKm1 = pathGrdK; pathGrdK = tempPtr;
- // 更新点索引
- idxOptVarStateKm1 = idxOptVarStateK;
- idxOptVarStateK = idxOptVarStateK + numState;
- idxOptVarControlKm1 = idxOptVarControlK;
- idxOptVarControlK = idxOptVarControlK + numControl;
- }
- if (dynRefK) { free(dynRefKm1); }
- if (dynGrdK) { free(dynGrdKm1); }
- if (pathRefK) { free(pathRefKm1); }
- if (pathGrdK) { free(pathGrdKm1); }
- if (dynRefK) { free(dynRefK); }
- if (dynGrdK) { free(dynGrdK); }
- if (pathRefK) { free(pathRefK); }
- if (pathGrdK) { free(pathGrdK); }
- // 更新段索引
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- free(dynGradInfo[0]);
- free(dynGradInfo);
- free(pathGradInfo[0]);
- free(pathGradInfo);
- socpInfo->cntNzAEqua = cntNzAEqua;
- socpInfo->cntRowAEqua = cntRowAEqua;
- socpInfo->cntNzAIneq = cntNzAIneq;
- socpInfo->cntRowAIneq = cntRowAIneq;
- socpInfo->cntNzASoc = cntNzASoc;
- socpInfo->cntRowASoc = cntRowASoc;
- socpInfo->cntSubASoc = cntSubASoc;
- }
- void socpInfoSparseLink(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_sparseInfo* sparseInfo, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup)
- {
- // 函数功能:获取衔接函数约束的二阶锥规划问题参数
- // 输入:
- // socpInfo:二阶锥信息结构体
- // refTraj:参考轨迹信息
- // consTypeInfo:约束类型信息
- // sparseInfo:稀疏(非零元位置)信息
- // setup:问题设置参数
- // 输出:
- // socpInfo:二阶锥信息结构体
- // 谢磊:2022 / 4 / 18编写
- // //基本参数
- int numPhase, * numPhaseOptVar, iphase, cnt, numNlEndp, * idxColAEqua, * ptrRowAEqua, * idxColACone, * ptrRowACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, cntSubASoc, nnzType1, i, p, j, i0, numType1, numType2, numType3, numType4, * idxType1, * idxType2, * idxType3, * idxType4, idxOptVarLeftPhaseStart, idxOptVarRightPhaseStart, numLeftMesh, numRightMesh, numLeftStatef, numLeftControl, numRightState, numRightControl, numLeftParameter, numNlLink, numSocLink, * dimSocLink, * leftStatefPtrRow, * leftTfPtrRow, * leftParameterPtrRow, * rightState0PtrRow, * rightT0PtrRow, * rightParameterPtrRow, * leftStatefIdxCol, * leftParameterIdxCol, * rightState0IdxCol, * rightParameterIdxCol, idxLinkGrdColLeftTf, idxLinkGrdColLeftParameter, idxLinkGrdColRightState0, idxLinkGrdColRightT0, idxLinkGrdColRightParameter, idxOptVarLeftStatef, idxOptVarLeftTf, idxOptVarLeftParameter, idxOptVarRightState0, idxOptVarRightT0, idxOptVarRightParameter,
- numLinkInput, numLinkOutput;
- double* valueAEqua, * bEqua, * valueACone, * bCone, temp, * nlLower, * nlUpper, * linkSocUpper, * refLeftStatef, refLeftTf, * refLeftParameter, * refRightState0, refRightT0, * refRightParameter, * linkRef, * linkGrd;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_bound* bounds;
- cmscp_consType** consType;
- cmscp_type* consTypePhase;
- cmscp_nzPos_link** linkNzPos;
- cmscp_derInfo** gradInfo;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- dimension = setup->dimension; // 问题维数信息
- mesh = setup->mesh; // 网格信息
- bounds = setup->bounds; // 上下边界信息
- numPhase = setup->numPhase; // 阶段数
- socpInfo = *socpInfo0;
- numPhaseOptVar = socpInfo->numPhaseOptVar; // 每段优化变量个数
- linkNzPos = sparseInfo->linkNzPos; // 稀疏信息
- consType = consTypeInfo->link; // 约束类型
- // 二阶锥规划参数
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- cntSubASoc = socpInfo->cntSubASoc;
- // 梯度信息
- gradInfo = (cmscp_derInfo**)malloc(numPhase * sizeof(cmscp_derInfo));
- gradInfoSparseLink(gradInfo, sparseInfo, refTraj, setup);
- //// 填充非零元素
- // 索引
- idxOptVarLeftPhaseStart = 0; // 各段起始索引
- idxOptVarRightPhaseStart = numPhaseOptVar[0]; // 各段起始索引
- for (iphase = 0; iphase < numPhase - 1; iphase++)
- {
-
- // 维数
- numLeftMesh = mesh[iphase]->num;
- numRightMesh = mesh[iphase + 1]->num;
- numLeftStatef = dimension->phase[iphase]->state;
- numLeftControl = dimension->phase[iphase]->control;
- numRightState = dimension->phase[iphase + 1]->state;
- numRightControl = dimension->phase[iphase + 1]->control;
- numLeftParameter = dimension->phase[iphase]->parameter;
- numNlLink = dimension->phase[iphase]->linkage->nl;
- numSocLink = dimension->phase[iphase]->linkage->soc->n;
- dimSocLink = dimension->phase[iphase]->linkage->soc->q;
- numLinkInput = numLeftStatef + numRightState + numRightState + numLeftParameter + 2;
- if (dimension->phase[iphase]->linkage->soc->n)
- {
- numLinkOutput = dimension->phase[iphase]->linkage->nl + cmscp_sum(dimension->phase[iphase]->linkage->soc->q, dimension->phase[iphase]->linkage->soc->n);
- }
- else
- {
- numLinkOutput = dimension->phase[iphase]->linkage->nl;
- }
- // 稀疏信息
- leftStatefPtrRow = linkNzPos[iphase]->leftStatef->ptrRow;
- leftTfPtrRow = linkNzPos[iphase]->leftTf->ptrRow;
- leftParameterPtrRow = linkNzPos[iphase]->leftParameter->ptrRow;
- rightState0PtrRow = linkNzPos[iphase]->rightState0->ptrRow;
- rightT0PtrRow = linkNzPos[iphase]->rightT0->ptrRow;
- rightParameterPtrRow = linkNzPos[iphase]->rightParameter->ptrRow;
- leftStatefIdxCol = linkNzPos[iphase]->leftStatef->idxCol;
- leftParameterIdxCol = linkNzPos[iphase]->leftParameter->idxCol;
- rightState0IdxCol = linkNzPos[iphase]->rightState0->idxCol;
- rightParameterIdxCol = linkNzPos[iphase]->rightParameter->idxCol;
- // 约束类型简写
- numType1 = consType[iphase]->type[0].num;
- numType2 = consType[iphase]->type[1].num;
- numType3 = consType[iphase]->type[2].num;
- numType4 = consType[iphase]->type[3].num;
- idxType1 = consType[iphase]->type[0].idx;
- idxType2 = consType[iphase]->type[1].idx;
- idxType3 = consType[iphase]->type[2].idx;
- idxType4 = consType[iphase]->type[3].idx;
- // 第1类约束的一半非零元个数
- nnzType1 = getType1LinkNnz(consType[iphase]->type[0], linkNzPos[iphase]);
- // 上下界
-
- nlLower = bounds->phase[iphase]->linkage->nl->lower;
- nlUpper = bounds->phase[iphase]->linkage->nl->upper;
- linkSocUpper = bounds->phase[iphase]->linkage->soc;
- // 参考轨迹
- refLeftStatef = refTraj[iphase]->state->v + (numLeftMesh - 1) * numLeftStatef;
- refLeftTf = refTraj[iphase]->finaltime;
- refLeftParameter = refTraj[iphase]->parameter;
- refRightState0 = refTraj[iphase + 1]->state->v;
- refRightT0 = refTraj[iphase + 1]->initialtime;
- refRightParameter = refTraj[iphase + 1]->parameter;
- // 梯度信息
- linkRef = gradInfo[iphase]->ref;
- linkGrd = gradInfo[iphase]->der;
- // 索引
- idxLinkGrdColLeftTf = numLeftStatef;
- idxLinkGrdColLeftParameter = idxLinkGrdColLeftTf + 1;
- idxLinkGrdColRightState0 = idxLinkGrdColLeftParameter + numLeftParameter;
- idxLinkGrdColRightT0 = idxLinkGrdColRightState0 + numRightState;
- idxLinkGrdColRightParameter = idxLinkGrdColRightT0 + 1;
- idxOptVarLeftStatef = idxOptVarLeftPhaseStart + (numLeftMesh - 1) * numLeftStatef;
- idxOptVarLeftTf = idxOptVarLeftPhaseStart + numLeftMesh * (numLeftStatef + numLeftControl) + 1;
- idxOptVarLeftParameter = idxOptVarLeftTf + 1;
- idxOptVarRightState0 = idxOptVarRightPhaseStart;
- idxOptVarRightT0 = idxOptVarRightPhaseStart + numRightMesh * (numRightState + numRightControl);
- idxOptVarRightParameter = idxOptVarRightT0 + 2;
- // 双不等式约束项,表达式:
- if (numType1)
- {
- // 上界:A * xfl + B * tfl + C * pl + D * x0r + E * t0r + F * pr <= b
- for (cnt = 0; cnt < numType1; cnt++) {
- i = idxType1[cnt];
- temp = -linkRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- ptrRowACone[cntRowAIneq + numType1] = cntNzAIneq + nnzType1;
- // 左侧xf前面系数
- for (p = leftStatefPtrRow[i]; p < leftStatefPtrRow[i + 1]; p++)
- {
- j = leftStatefIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarLeftStatef + j;
- valueACone[cntNzAIneq] = linkGrd[i + j* numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refLeftStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 左侧tf前面系数
- if (leftTfPtrRow[i] + 1 == leftTfPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarLeftTf;
- valueACone[cntNzAIneq] = linkGrd[i + idxLinkGrdColLeftTf * numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refLeftTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 左侧p前面系数
- for (p = leftParameterPtrRow[i]; p < leftParameterPtrRow[i + 1]; p++)
- {
- j = leftParameterIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarLeftParameter + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColLeftParameter + j) * numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refLeftParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧x0前面系数
- for (p = rightState0PtrRow[i]; p < rightState0PtrRow[i + 1]; p++)
- {
- j = rightState0IdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarRightState0 + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColRightState0 + j) * numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refRightState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧t0前面系数
- if (rightT0PtrRow[i] + 1 == rightT0PtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarRightT0;
- valueACone[cntNzAIneq] = linkGrd[i + idxLinkGrdColRightT0 * numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refRightT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧p前面系数
- for (p = rightParameterPtrRow[i]; p < rightParameterPtrRow[i + 1]; p++) {
- j = rightParameterIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarRightParameter + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColRightParameter + j) * numLinkOutput];
- idxColACone[cntNzAIneq + nnzType1] = idxColACone[cntNzAIneq];
- valueACone[cntNzAIneq + nnzType1] = -valueACone[cntNzAIneq];
- temp = temp + valueACone[cntNzAIneq] * refRightParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp + nlUpper[i];
- bCone[cntRowAIneq] = -temp - nlLower[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- // 等式约束项,表达式:
- // A* xfl + B * tfl + C * pl + D * x0r + E * t0r + F * pr == b
- if (numType2)
- {
- for (cnt = 0; cnt < numType2; cnt++) {
- i = idxType2[cnt];
- temp = -linkRef[i];
- // 每行首个非零元的位置
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- // 左侧xf前面系数
- for (p = leftStatefPtrRow[i]; p < leftStatefPtrRow[i + 1]; p++)
- {
- j = leftStatefIdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarLeftStatef + j;
- valueAEqua[cntNzAEqua] = linkGrd[i + j * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refLeftStatef[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // 左侧tf前面系数
- if (leftTfPtrRow[i] + 1 == leftTfPtrRow[i + 1])
- {
- idxColAEqua[cntNzAEqua] = idxOptVarLeftTf;
- valueAEqua[cntNzAEqua] = linkGrd[i + (idxLinkGrdColLeftTf) * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refLeftTf;
- cntNzAEqua = cntNzAEqua + 1;
- }
- // 左侧p前面系数
- for (p = leftParameterPtrRow[i]; p < leftParameterPtrRow[i + 1]; p++)
- {
- j = leftParameterIdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarLeftParameter + j;
- valueAEqua[cntNzAEqua] = linkGrd[i + (idxLinkGrdColLeftParameter + j) * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refLeftParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // 右侧x0前面系数
- for (p = rightState0PtrRow[i]; p < rightState0PtrRow[i + 1]; p++) {
- j = rightState0IdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarRightState0 + j;
- valueAEqua[cntNzAEqua] = linkGrd[i + (idxLinkGrdColRightState0 + j) * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refRightState0[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- // 右侧t0前面系数
- if (rightT0PtrRow[i] + 1 == rightT0PtrRow[i + 1]) {
- idxColAEqua[cntNzAEqua] = idxOptVarRightT0;
- valueAEqua[cntNzAEqua] = linkGrd[i + idxLinkGrdColRightT0 * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refRightT0;
- cntNzAEqua = cntNzAEqua + 1;
- }
- // 右侧p前面系数
- for (p = rightParameterPtrRow[i]; p < rightParameterPtrRow[i + 1]; p++)
- {
- j = rightParameterIdxCol[p]; // 列
- idxColAEqua[cntNzAEqua] = idxOptVarRightParameter + j;
- valueAEqua[cntNzAEqua] = linkGrd[i + (idxLinkGrdColRightParameter + j) * numLinkOutput];
- temp = temp + valueAEqua[cntNzAEqua] * refRightParameter[j];
- cntNzAEqua = cntNzAEqua + 1;
- }
- bEqua[cntRowAEqua] = temp + nlUpper[i];
- cntRowAEqua = cntRowAEqua + 1;
- }
- }
- // 上不等式约束项,表达式:
- // A* xfl + B * tfl + C * pl + D * x0r + E * t0r + F * pr <= b
- if (numType3)
- {
- for (cnt = 0; cnt < numType3; cnt++) {
- i = idxType3[cnt];
- temp = -linkRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // 左侧xf前面系数
- for (p = leftStatefPtrRow[i]; p = leftStatefPtrRow[i + 1]; p++) {
- j = leftStatefIdxCol[p]; // 列
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarLeftStatef + j;
- valueACone[cntNzAIneq] = linkGrd[i + j * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftStatef[j];
- }
- // 左侧tf前面系数
- if (leftTfPtrRow[i] + 1 == leftTfPtrRow[i + 1])
- {
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarLeftTf;
- valueACone[cntNzAIneq] = linkGrd[i + idxLinkGrdColLeftTf * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftTf;
- }
- // 左侧p前面系数
- for (p = leftParameterPtrRow[i]; p < leftParameterPtrRow[i + 1]; p++) {
- j = leftParameterIdxCol[p]; // 列
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarLeftParameter + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColLeftParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftParameter[j];
- }
- // 右侧x0前面系数
- for (p = rightState0PtrRow[i]; p < rightState0PtrRow[i + 1]; p++) {
- j = rightState0IdxCol[p]; // 列
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarRightState0 + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColRightState0 + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightState0[j];
- }
- // 右侧t0前面系数
- if (rightT0PtrRow[i] + 1 == rightT0PtrRow[i + 1])
- {
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarRightT0 + 1;
- valueACone[cntNzAIneq] = linkGrd[i + idxLinkGrdColRightT0 * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightT0;
- }
- // 右侧p前面系数
- for (p = rightParameterPtrRow[i]; p < rightParameterPtrRow[i + 1]; p++)
- {
- j = rightParameterIdxCol[p]; // 列
- cntNzAIneq = cntNzAIneq + 1;
- idxColACone[cntNzAIneq] = idxOptVarRightParameter + j;
- valueACone[cntNzAIneq] = linkGrd[i + (idxLinkGrdColRightParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightParameter[j];
- }
- bCone[cntRowAIneq] = temp + nlUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- // 下不等式约束项,表达式:
- // -A * xfl - B * tfl - C * pl - D * x0r - E * t0r - F * pr <= -b
- if (numType4)
- {
- for (cnt = 0; cnt < numType4; cnt++)
- {
- i = idxType4[cnt];
- temp = linkRef[i];
- // 每行首个非零元的位置
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- // 左侧xf前面系数
- for (p = leftStatefPtrRow[i]; p < leftStatefPtrRow[i + 1]; p++)
- {
- j = leftStatefIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarLeftStatef + j;
- valueACone[cntNzAIneq] = -linkGrd[i + j * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftStatef[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 左侧tf前面系数
- if (leftTfPtrRow[i] + 1 == leftTfPtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarLeftTf;
- valueACone[cntNzAIneq] = -linkGrd[i + idxLinkGrdColLeftTf * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftTf;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 左侧p前面系数
- for (p = leftParameterPtrRow[i]; p < leftParameterPtrRow[i + 1]; p++)
- {
- j = leftParameterIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarLeftParameter + j;
- valueACone[cntNzAIneq] = -linkGrd[i + (idxLinkGrdColLeftParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refLeftParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧x0前面系数
- for (p = rightState0PtrRow[i]; p < rightState0PtrRow[i + 1]; p++)
- {
- j = rightState0IdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarRightState0 + j;
- valueACone[cntNzAIneq] = -linkGrd[i + (idxLinkGrdColRightState0 + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightState0[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧t0前面系数
- if (rightT0PtrRow[i] + 1 == rightT0PtrRow[i + 1])
- {
- idxColACone[cntNzAIneq] = idxOptVarRightT0;
- valueACone[cntNzAIneq] = -linkGrd[i + idxLinkGrdColRightT0 * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightT0;
- cntNzAIneq = cntNzAIneq + 1;
- }
- // 右侧p前面系数
- for (p = rightParameterPtrRow[i]; p < rightParameterPtrRow[i + 1]; p++)
- {
- j = rightParameterIdxCol[p]; // 列
- idxColACone[cntNzAIneq] = idxOptVarRightParameter + j;
- valueACone[cntNzAIneq] = -linkGrd[i + (idxLinkGrdColRightParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzAIneq] * refRightParameter[j];
- cntNzAIneq = cntNzAIneq + 1;
- }
- bCone[cntRowAIneq] = temp - nlLower[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- }
- if (numSocLink)
- {
- for (i = 0; i < cmscp_sum(dimSocLink, numSocLink); i++)
- {
- cntRowASoc = cntRowASoc + 1;
- i0 = i + numNlLink;
- temp = -linkRef[i0];
- ptrRowACone[cntRowASoc] = cntNzASoc;
- // 左侧xf前面系数
- for (p = leftStatefPtrRow[i0]; p < leftStatefPtrRow[i0 + 1]; p++)
- {
- j = leftStatefIdxCol[p]; // 列
- idxColACone[cntNzASoc] = idxOptVarLeftStatef + j;
- valueACone[cntNzASoc] = linkGrd[i0 + j * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refLeftStatef[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // 左侧tf前面系数
- if (leftTfPtrRow[i0] + 1 == leftTfPtrRow[i0 + 1])
- {
- idxColACone[cntNzASoc] = idxOptVarLeftTf;
- valueACone[cntNzASoc] = linkGrd[i0 + idxLinkGrdColLeftTf * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refLeftTf;
- cntNzASoc = cntNzASoc + 1;
- }
- // 左侧p前面系数
- for (p = leftParameterPtrRow[i0]; p < leftParameterPtrRow[i0 + 1]; p++)
- {
- j = leftParameterIdxCol[p]; // 列
- idxColACone[cntNzASoc] = idxOptVarLeftParameter + j;
- valueACone[cntNzASoc] = linkGrd[i0 + (idxLinkGrdColLeftParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refLeftParameter[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // 右侧x0前面系数
- for (p = rightState0PtrRow[i0]; p < rightState0PtrRow[i0 + 1]; p++)
- {
- j = rightState0IdxCol[p]; // 列
- idxColACone[cntNzASoc] = idxOptVarRightState0 + j;
- valueACone[cntNzASoc] = linkGrd[i0 + (idxLinkGrdColRightState0 + j) * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refRightState0[j];
- cntNzASoc = cntNzASoc + 1;
- }
- // 右侧t0前面系数
- if (rightT0PtrRow[i0] + 1 == rightT0PtrRow[i0 + 1])
- {
- idxColACone[cntNzASoc] = idxOptVarRightT0 + 1;
- valueACone[cntNzASoc] = linkGrd[i0 + idxLinkGrdColRightT0 * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refRightT0;
- cntNzASoc = cntNzASoc + 1;
- }
- // 右侧p前面系数
- for (p = rightParameterPtrRow[i0]; p < rightParameterPtrRow[i0 + 1]; p++)
- {
- j = rightParameterIdxCol[p]; // 列
- idxColACone[cntNzASoc] = idxOptVarRightParameter + j;
- valueACone[cntNzASoc] = linkGrd[i0 + (idxLinkGrdColRightParameter + j) * numLinkOutput];
- temp = temp + valueACone[cntNzASoc] * refRightParameter[j];
- cntNzASoc = cntNzASoc + 1;
- }
- bCone[cntRowASoc] = temp + linkSocUpper[i];
- cntRowAIneq = cntRowAIneq + 1;
- }
- for (i = 0; i < numSocLink; i++)
- {
- cntSubASoc = cntSubASoc + 1;
- dimCone->q[cntSubASoc] = dimSocLink[i];
- }
- }
- idxOptVarLeftPhaseStart = idxOptVarRightPhaseStart;
- idxOptVarRightPhaseStart = idxOptVarRightPhaseStart + numPhaseOptVar[iphase + 1];
- }
- for ( iphase = 0; iphase < numPhase-1; iphase++)
- {
- free(gradInfo[iphase]->ref);
- free(gradInfo[iphase]->der);
- free(gradInfo[iphase]);
- }
- free(gradInfo);
- socpInfo->cntNzAEqua = cntNzAEqua;
- socpInfo->cntRowAEqua = cntRowAEqua;
- socpInfo->cntNzAIneq = cntNzAIneq;
- socpInfo->cntRowAIneq = cntRowAIneq;
- socpInfo->cntNzASoc = cntNzASoc;
- socpInfo->cntRowASoc = cntRowASoc;
- socpInfo->cntSubASoc = cntSubASoc;
- }
- void socpInfoBound(cmscp_socpInfo** socpInfo0, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup) {
- // socpInfo: 已有的二阶锥信息
- // refTraj : 参考轨迹信息
- // gradInfo : 梯度信息
- // setup : 设置信息
- int numPhase, * numPhaseOptVar, iphase, numMesh, numState, numControl, cnt, * idxColAEqua, * ptrRowAEqua, * idxColACone, * ptrRowACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, idxOptVarPhaseStart, idxOptVarState0, idxOptVarStatef, idxOptVarT0, idxOptVarTf, idxOptVarParameter, i, idxOptVarStateK, idxOptVarControlK, k, numOptVar, numIniTimeType1, numIniTimeType2, numIniTimeType3, numIniTimeType4, * idxIniTimeType1, * idxIniTimeType2, * idxIniTimeType3, * idxIniTimeType4, numFinTimeType1, numFinTimeType2, numFinTimeType3, numFinTimeType4, * idxFinTimeType1, * idxFinTimeType2, * idxFinTimeType3, * idxFinTimeType4, numIniStateType1, numIniStateType2, numIniStateType3, numIniStateType4, * idxIniStateType1, * idxIniStateType2, * idxIniStateType3, * idxIniStateType4, numFinStateType1, numFinStateType2, numFinStateType3, numFinStateType4, * idxFinStateType1, * idxFinStateType2, * idxFinStateType3, * idxFinStateType4, numStateType1, numStateType2, numStateType3, numStateType4, * idxStateType1, * idxStateType2, * idxStateType3, * idxStateType4, numControlType1, numControlType2, numControlType3, numControlType4, * idxControlType1, * idxControlType2, * idxControlType3, * idxControlType4, numParameterType1, numParameterType2, numParameterType3, numParameterType4, * idxParameterType1, * idxParameterType2, * idxParameterType3, * idxParameterType4;
- double* valueAEqua, * bEqua, * valueACone, * bCone, iniTimeLower, iniTimeUpper, finTimeLower, finTimeUpper, * iniStateLower, * iniStateUpper, * finStateLower, * finStateUpper, * stateLower, * stateUpper, * controlLower, * controlUpper, * parameterLower, * parameterUpper;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_bound* bounds;
- cmscp_consType_bound* consType;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- cmscp_consType** iniTimeType, ** finTimeType, ** iniStateType, ** finStateType, ** stateType, ** controlType, ** parameterType;
- // 1 基本参数
- dimension = setup->dimension;// 问题维数信息
- mesh = setup->mesh;// 网格信息
- bounds = setup->bounds;// 上下边界信息
- numPhase = setup->numPhase;// 阶段数
- socpInfo = *socpInfo0;
- numOptVar = socpInfo->numOptVar;// 总优化变量个数
- numPhaseOptVar = socpInfo->numPhaseOptVar;// 每段优化变量个数
- consType = consTypeInfo->bound;// 约束类型
- iniTimeType = consType->iniTime;
- finTimeType = consType->finTime;
- iniStateType = consType->iniState;
- finStateType = consType->finState;
- stateType = consType->state;
- controlType = consType->control;
- parameterType = consType->parameter;
- // 二阶锥规划参数
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- // 2 填充非零元素
- // 各段起始变量索引
- idxOptVarPhaseStart = 0;
- for (iphase = 0; iphase < numPhase; iphase++) {
- // 维数
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- // 初始时刻约束类型简写
- numIniTimeType1 = iniTimeType[iphase]->type[0].num;
- numIniTimeType2 = iniTimeType[iphase]->type[1].num;
- numIniTimeType3 = iniTimeType[iphase]->type[2].num;
- numIniTimeType4 = iniTimeType[iphase]->type[3].num;
- idxIniTimeType1 = iniTimeType[iphase]->type[0].idx;
- idxIniTimeType2 = iniTimeType[iphase]->type[1].idx;
- idxIniTimeType3 = iniTimeType[iphase]->type[2].idx;
- idxIniTimeType4 = iniTimeType[iphase]->type[3].idx;
- // 终端时刻约束类型简写
- numFinTimeType1 = finTimeType[iphase]->type[0].num;
- numFinTimeType2 = finTimeType[iphase]->type[1].num;
- numFinTimeType3 = finTimeType[iphase]->type[2].num;
- numFinTimeType4 = finTimeType[iphase]->type[3].num;
- idxFinTimeType1 = finTimeType[iphase]->type[0].idx;
- idxFinTimeType2 = finTimeType[iphase]->type[1].idx;
- idxFinTimeType3 = finTimeType[iphase]->type[2].idx;
- idxFinTimeType4 = finTimeType[iphase]->type[3].idx;
- // 初始状态约束类型简写
- numIniStateType1 = iniStateType[iphase]->type[0].num;
- numIniStateType2 = iniStateType[iphase]->type[1].num;
- numIniStateType3 = iniStateType[iphase]->type[2].num;
- numIniStateType4 = iniStateType[iphase]->type[3].num;
- idxIniStateType1 = iniStateType[iphase]->type[0].idx;
- idxIniStateType2 = iniStateType[iphase]->type[1].idx;
- idxIniStateType3 = iniStateType[iphase]->type[2].idx;
- idxIniStateType4 = iniStateType[iphase]->type[3].idx;
- // 终端状态约束类型简写
- numFinStateType1 = finStateType[iphase]->type[0].num;
- numFinStateType2 = finStateType[iphase]->type[1].num;
- numFinStateType3 = finStateType[iphase]->type[2].num;
- numFinStateType4 = finStateType[iphase]->type[3].num;
- idxFinStateType1 = finStateType[iphase]->type[0].idx;
- idxFinStateType2 = finStateType[iphase]->type[1].idx;
- idxFinStateType3 = finStateType[iphase]->type[2].idx;
- idxFinStateType4 = finStateType[iphase]->type[3].idx;
- // 状态约束类型简写
- numStateType1 = stateType[iphase]->type[0].num;
- numStateType2 = stateType[iphase]->type[1].num;
- numStateType3 = stateType[iphase]->type[2].num;
- numStateType4 = stateType[iphase]->type[3].num;
- idxStateType1 = stateType[iphase]->type[0].idx;
- idxStateType2 = stateType[iphase]->type[1].idx;
- idxStateType3 = stateType[iphase]->type[2].idx;
- idxStateType4 = stateType[iphase]->type[3].idx;
- // 控制约束类型简写
- numControlType1 = controlType[iphase]->type[0].num;
- numControlType2 = controlType[iphase]->type[1].num;
- numControlType3 = controlType[iphase]->type[2].num;
- numControlType4 = controlType[iphase]->type[3].num;
- idxControlType1 = controlType[iphase]->type[0].idx;
- idxControlType2 = controlType[iphase]->type[1].idx;
- idxControlType3 = controlType[iphase]->type[2].idx;
- idxControlType4 = controlType[iphase]->type[3].idx;
- // 静态变量约束类型简写
- numParameterType1 = parameterType[iphase]->type[0].num;
- numParameterType2 = parameterType[iphase]->type[1].num;
- numParameterType3 = parameterType[iphase]->type[2].num;
- numParameterType4 = parameterType[iphase]->type[3].num;
- idxParameterType1 = parameterType[iphase]->type[0].idx;
- idxParameterType2 = parameterType[iphase]->type[1].idx;
- idxParameterType3 = parameterType[iphase]->type[2].idx;
- idxParameterType4 = parameterType[iphase]->type[3].idx;
- // 初始和终端时刻约束边界简写
- iniTimeLower = bounds->phase[iphase]->initialtime->lower;
- iniTimeUpper = bounds->phase[iphase]->initialtime->upper;
- finTimeLower = bounds->phase[iphase]->finaltime->lower;
- finTimeUpper = bounds->phase[iphase]->finaltime->upper;
- // 初始、终端和过程状态约束边界简写
- iniStateLower = bounds->phase[iphase]->initialstate->lower;
- iniStateUpper = bounds->phase[iphase]->initialstate->upper;
- finStateLower = bounds->phase[iphase]->finalstate->lower;
- finStateUpper = bounds->phase[iphase]->finalstate->upper;
- stateLower = bounds->phase[iphase]->state->lower;
- stateUpper = bounds->phase[iphase]->state->upper;
- // 控制约束边界简写
- controlLower = bounds->phase[iphase]->control->lower;
- controlUpper = bounds->phase[iphase]->control->upper;
- // 静态变量约束边界简写
- parameterLower = bounds->phase[iphase]->parameter->lower;
- parameterUpper = bounds->phase[iphase]->parameter->upper;
- // 各个变量索引
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarTf + 1;
- // 初始状态变量
- idxOptVarState0 = idxOptVarPhaseStart;
- if (numIniStateType1) {
- for (cnt = 0; cnt < numIniStateType1; cnt++) {
- i = idxIniStateType1[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = iniStateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarState0 + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -xK <= -bLower
- bCone[cntRowAIneq] = -iniStateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarState0 + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numIniStateType2) {
- for (cnt = 0; cnt < numIniStateType2; cnt++) {
- i = idxIniStateType2[cnt];
- // xK == bUpper
- bEqua[cntRowAEqua] = iniStateUpper[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarState0 + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- }
- if (numIniStateType3) {
- for (cnt = 0; cnt < numIniStateType3; cnt++) {
- i = idxIniStateType3[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = iniStateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarState0 + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numIniStateType4) {
- for (cnt = 0; cnt < numIniStateType4; cnt++) {
- i = idxIniStateType4[cnt];
- // -xK <= -bLower
- bCone[cntRowAIneq] = -iniStateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarState0 + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- // 过程状态
- idxOptVarStateK = idxOptVarPhaseStart + numState;
- if (numStateType1 != 0) {
- for (k = 1; k < numMesh - 1; k++) {
- for (cnt = 0; cnt < numStateType1; cnt++) {
- i = idxStateType1[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = stateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStateK + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -xK <= -bLower
- bCone[cntRowAIneq] = -stateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStateK + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarStateK = idxOptVarStateK + numState;
- }
- }
- idxOptVarStateK = idxOptVarPhaseStart + numState;
- if (numStateType2 != 0) {
- for (k = 1; k < numMesh - 1; k++) {
- for (cnt = 0; cnt < numStateType2; cnt++) {
- i = idxStateType2[cnt];
- // xK == bUpper
- bEqua[cntRowAEqua] = stateUpper[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarStateK + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxOptVarStateK = idxOptVarStateK + numState;
- }
- }
- idxOptVarStateK = idxOptVarPhaseStart + numState;
- if (numStateType3 != 0) {
- for (k = 1; k < numMesh - 1; k++) {
- for (cnt = 0; cnt < numStateType3; cnt++) {
- i = idxStateType3[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = stateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStateK + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarStateK = idxOptVarStateK + numState;
- }
- }
- idxOptVarStateK = idxOptVarPhaseStart + numState;
- if (numStateType4 != 0) {
- for (k = 1; k < numMesh - 1; k++) {
- for (cnt = 0; cnt < numStateType4; cnt++) {
- i = idxStateType4[cnt];
- // -xK <= -bLower
- bCone[cntRowAIneq] = -stateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStateK + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarStateK = idxOptVarStateK + numState;
- }
- }
- // 终端状态
- idxOptVarStatef = idxOptVarPhaseStart + (numMesh - 1) * numState;
- if (numFinStateType1 != 0) {
- for (cnt = 0; cnt < numFinStateType1; cnt++) {
- i = idxFinStateType1[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = finStateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStatef + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -xK <= -bLower
- bCone[cntRowAIneq] = -finStateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStatef + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numFinStateType2 != 0) {
- for (cnt = 0; cnt < numFinStateType2; cnt++) {
- i = idxFinStateType2[cnt];
- // xK == bUpper
- bEqua[cntRowAEqua] = finStateUpper[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarStatef + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- }
- if (numFinStateType3 != 0) {
- for (cnt = 0; cnt < numFinStateType3; cnt++) {
- i = idxFinStateType3[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = finStateUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStatef + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numFinStateType4 != 0) {
- for (cnt = 0; cnt < numFinStateType4; cnt++) {
- i = idxFinStateType4[cnt];
- // -xK <= -bLower
- bCone[cntRowAIneq] = -finStateLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarStatef + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- // 控制变量
- idxOptVarControlK = idxOptVarPhaseStart + numMesh * numState;
- if (numControlType1 != 0) {
- for (k = 0; k < numMesh; k++) {
- for (cnt = 0; cnt < numControlType1; cnt++) {
- i = idxControlType1[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = controlUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarControlK + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -xK <= -bLower
- bCone[cntRowAIneq] = -controlLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarControlK + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarControlK = idxOptVarControlK + numControl;
- }
- }
- idxOptVarControlK = idxOptVarPhaseStart + numMesh * numState;
- if (numControlType2 != 0) {
- for (k = 0; k < numMesh; k++) {
- for (cnt = 0; cnt < numControlType2; cnt++) {
- i = idxControlType2[cnt];
- // xK == bUpper
- bEqua[cntRowAEqua] = controlUpper[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarControlK + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- idxOptVarControlK = idxOptVarControlK + numControl;
- }
- }
- idxOptVarControlK = idxOptVarPhaseStart + numMesh * numState;
- if (numControlType3 != 0) {
- for (k = 0; k < numMesh; k++) {
- for (cnt = 0; cnt < numControlType3; cnt++) {
- i = idxControlType3[cnt];
- // xK <= bUpper
- bCone[cntRowAIneq] = controlUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarControlK + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarControlK = idxOptVarControlK + numControl;
- }
- }
- idxOptVarControlK = idxOptVarPhaseStart + numMesh * numState;
- if (numControlType4 != 0) {
- for (k = 0; k < numMesh; k++) {
- for (cnt = 0; cnt < numControlType4; cnt++) {
- i = idxControlType4[cnt];
- // -xK <= -bLower
- bCone[cntRowAIneq] = -controlLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarControlK + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- idxOptVarControlK = idxOptVarControlK + numControl;
- }
- }
- // 初始时刻
- if (numIniTimeType1 != 0) {
- for (cnt = 0; cnt < numIniTimeType1; cnt++) {
- i = idxIniTimeType1[cnt];
- // t0 <= bUpper
- bCone[cntRowAIneq] = iniTimeUpper;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarT0 + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -t0 <= -bLower
- bCone[cntRowAIneq] = -iniTimeLower;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarT0 + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numIniTimeType2 != 0) {
- for (cnt = 0; cnt < numIniTimeType2; cnt++) {
- i = idxIniTimeType2[cnt];
- // t0 == bUpper
- bEqua[cntRowAEqua] = iniTimeUpper;
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarT0 + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- }
- if (numIniTimeType3 != 0) {
- for (cnt = 0; cnt < numIniTimeType3; cnt++) {
- i = idxIniTimeType3[cnt];
- // t0 <= bUpper
- bCone[cntRowAIneq] = iniTimeUpper;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarT0 + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numIniTimeType4 != 0) {
- for (cnt = 0; cnt < numIniTimeType4; cnt++) {
- i = idxIniTimeType4[cnt];
- // -t0 <= -bLower
- bCone[cntRowAIneq] = -iniTimeLower;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarT0 + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- // 终端时刻
- if (numFinTimeType1 != 0) {
- for (cnt = 0; cnt < numFinTimeType1; cnt++) {
- i = idxFinTimeType1[cnt];
- // tf <= bUpper
- bCone[cntRowAIneq] = finTimeUpper;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarTf + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -tf <= -bLower
- bCone[cntRowAIneq] = -finTimeLower;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarTf + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numFinTimeType2 != 0) {
- for (cnt = 0; cnt < numFinTimeType2; cnt++) {
- i = idxFinTimeType2[cnt];
- // tf == bUpper
- bEqua[cntRowAEqua] = finTimeUpper;
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarTf + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- }
- if (numFinTimeType3 != 0) {
- for (cnt = 0; cnt < numFinTimeType3; cnt++) {
- i = idxFinTimeType3[cnt];
- // tf <= bUpper
- bCone[cntRowAIneq] = finTimeUpper;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarTf + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numFinTimeType4 != 0) {
- for (cnt = 0; cnt < numFinTimeType4; cnt++) {
- i = idxFinTimeType4[cnt];
- // -tf <= -bLower
- bCone[cntRowAIneq] = -finTimeLower;
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarTf + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- // 静态参数变量
- if (numParameterType1 != 0) {
- for (cnt = 0; cnt < numParameterType1; cnt++) {
- i = idxParameterType1[cnt];
- // p <= bUpper
- bCone[cntRowAIneq] = parameterUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarParameter + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- // -p <= -bLower
- bCone[cntRowAIneq] = -parameterLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarParameter + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numParameterType2 != 0) {
- for (cnt = 0; cnt < numParameterType2; cnt++) {
- i = idxParameterType2[cnt];
- // p == bUpper
- bEqua[cntRowAEqua] = parameterUpper[i];
- ptrRowAEqua[cntRowAEqua] = cntNzAEqua;
- idxColAEqua[cntNzAEqua] = idxOptVarParameter + i;
- valueAEqua[cntNzAEqua] = 1;
- cntRowAEqua = cntRowAEqua + 1;
- cntNzAEqua = cntNzAEqua + 1;
- }
- }
- if (numParameterType3 != 0) {
- for (cnt = 0; cnt < numParameterType3; cnt++) {
- i = idxParameterType3[cnt];
- // p <= bUpper
- bCone[cntRowAIneq] = parameterUpper[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarParameter + i;
- valueACone[cntNzAIneq] = 1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- if (numParameterType4 != 0) {
- for (cnt = 0; cnt < numParameterType4; cnt++) {
- i = idxParameterType4[cnt];
- // -p <= -bLower
- bCone[cntRowAIneq] = -parameterLower[i];
- ptrRowACone[cntRowAIneq] = cntNzAIneq;
- idxColACone[cntNzAIneq] = idxOptVarParameter + i;
- valueACone[cntNzAIneq] = -1.0;
- cntRowAIneq = cntRowAIneq + 1;
- cntNzAIneq = cntNzAIneq + 1;
- }
- }
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- socpInfo->cntNzAEqua = cntNzAEqua;
- socpInfo->cntRowAEqua = cntRowAEqua;
- socpInfo->cntNzAIneq = cntNzAIneq;
- socpInfo->cntRowAIneq = cntRowAIneq;
- socpInfo->cntNzASoc = cntNzASoc;
- socpInfo->cntRowASoc = cntRowASoc;
- }
- void socpInfoTrust(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_setup* setup)
- {
- // socpInfo: 已有的二阶锥信息
- // refTraj : 参考轨迹信息
- // gradInfo : 梯度信息
- // setup : 设置信息
- //
- int numPhase, * numPhaseOptVar, iphase, numMesh, numState, numControl, * idxColAEqua, * ptrRowAEqua, * idxColACone, * ptrRowACone, cntNzAEqua, cntRowAEqua, cntNzAIneq, cntRowAIneq, cntNzASoc, cntRowASoc, cntSubASoc, idxOptVarPhaseStart, idxOptVarT0, idxOptVarTf, idxOptVarParameter, p, j, numParameter, k, idxStateK, idxControlK, idxOptVarTrust;
- double* valueAEqua, * bEqua, * valueACone, * bCone, refT0, refTf, * refParameter, * refState, * refControl, a1, a2, weightIniTime, weightFinTime, * weightState, * weightControl, * weightParameter;
- cmscp_dimension* dimension;
- cmscp_mesh** mesh;
- cmscp_trust* trust;
- cmscp_socpInfo* socpInfo;
- cmscp_dimCone* dimCone;
- // //基本参数
- dimension = setup->dimension;// 问题维数信息
- mesh = setup->mesh;// 网格信息
- numPhase = setup->numPhase;// 阶段数
- socpInfo = *socpInfo0;
- numPhaseOptVar = socpInfo->numPhaseOptVar;// 每段优化变量个数
- trust = setup->trust;// 信赖域信息
- // 二阶锥规划问题参数
- // 二阶锥规划参数
- ptrRowAEqua = socpInfo->AEqua->ptrRow;
- idxColAEqua = socpInfo->AEqua->idxCol;
- valueAEqua = socpInfo->AEqua->v;
- bEqua = socpInfo->bEqua;
- ptrRowACone = socpInfo->ACone->ptrRow;
- idxColACone = socpInfo->ACone->idxCol;
- valueACone = socpInfo->ACone->v;
- bCone = socpInfo->bCone;
- dimCone = socpInfo->dimCone;
- // 计数变量
- cntNzAEqua = socpInfo->cntNzAEqua;
- cntRowAEqua = socpInfo->cntRowAEqua;
- cntNzAIneq = socpInfo->cntNzAIneq;
- cntRowAIneq = socpInfo->cntRowAIneq;
- cntNzASoc = socpInfo->cntNzASoc;
- cntRowASoc = socpInfo->cntRowASoc;
- cntSubASoc = socpInfo->cntSubASoc;
- //// 填充信赖域约束非零元
- // 硬信赖域:
- // [0] <= _K[delta]
- // [x][x_r]
- // 或者软信赖域:
- // [-delta] <= _K[0]
- // [x][x_r]
- // 每段起始点索引
- idxOptVarPhaseStart = 0;
- for (iphase = 0; iphase < numPhase; iphase++) {
- numMesh = mesh[iphase]->num;
- numState = dimension->phase[iphase]->state;
- numControl = dimension->phase[iphase]->control;
- numParameter = dimension->phase[iphase]->parameter;
- refState = refTraj[iphase]->state->v;
- refControl = refTraj[iphase]->control->v;
- refT0 = refTraj[iphase]->initialtime;
- refTf = refTraj[iphase]->finaltime;
- refParameter = refTraj[iphase]->parameter;
- if (strcmp(trust->model, "hard")==0) {
- a1 = 0;
- a2 = trust->phase[iphase]->hardRadius;
- }
- else if (strcmp(trust->model, "soft")==0) {
- a1 = -1;
- a2 = 0;
- }
- weightIniTime = trust->phase[iphase]->weight->initialtime;
- weightFinTime = trust->phase[iphase]->weight->finaltime;
- weightState = trust->phase[iphase]->weight->state;
- weightControl = trust->phase[iphase]->weight->control;
- weightParameter = trust->phase[iphase]->weight->parameter;
- // 索引
-
- idxStateK = idxOptVarPhaseStart;
- idxControlK = idxOptVarPhaseStart + numMesh * numState;
- idxOptVarT0 = idxOptVarPhaseStart + numMesh * (numState + numControl);
- idxOptVarTf = idxOptVarT0 + 1;
- idxOptVarParameter = idxOptVarTf + 1;
- idxOptVarTrust = idxOptVarParameter + numParameter;
- for (k = 0; k < numMesh; k++) {
- // 状态变量
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTrust;
- valueACone[cntNzASoc] = a1;
- bCone[cntRowASoc] = a2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- for (j = 0; j < numState; j++) {
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxStateK + j;
- valueACone[cntNzASoc] = weightState[j];
- bCone[cntRowASoc] = weightState[j] * refState[j + k * numState];
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- }
- dimCone->q[cntSubASoc] = numState + 1;
- cntSubASoc = cntSubASoc + 1;
- // 控制变量
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTrust;
- valueACone[cntNzASoc] = a1;
- bCone[cntRowASoc] = a2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- for (j = 0; j < numControl; j++) {
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxControlK + j;
- valueACone[cntNzASoc] = weightControl[j];
- bCone[cntRowASoc] = weightControl[j] * refControl[j + k* numControl];
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- }
- dimCone->q[cntSubASoc] = numControl + 1;
- cntSubASoc = cntSubASoc + 1;
- idxStateK = idxStateK + numState;
- idxControlK = idxControlK + numControl;
- }
- // 初始时刻
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTrust;
- valueACone[cntNzASoc] = a1;
- bCone[cntRowASoc] = a2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarT0;
- valueACone[cntNzASoc] = weightIniTime;
- bCone[cntRowASoc] = weightIniTime * refT0;
- dimCone->q[cntSubASoc] = 2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- cntSubASoc = cntSubASoc + 1;
- // 终端时刻
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTrust;
- valueACone[cntNzASoc] = a1;
- bCone[cntRowASoc] = a2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTf;
- valueACone[cntNzASoc] = weightFinTime;
- bCone[cntRowASoc] = weightFinTime * refTf;
- dimCone->q[cntSubASoc] = 2;
- cntSubASoc = cntSubASoc + 1;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- // 参数变量
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarTrust;
- valueACone[cntNzASoc] = a1;
- bCone[cntRowASoc] = a2;
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- for (j = 0; j < numParameter; j++) {
- ptrRowACone[cntRowASoc] = cntNzASoc;
- idxColACone[cntNzASoc] = idxOptVarParameter + j;
- valueACone[cntNzASoc] = weightParameter[j];
- bCone[cntRowASoc] = weightParameter[j] * refParameter[j];
- cntRowASoc = cntRowASoc + 1;
- cntNzASoc = cntNzASoc + 1;
- }
- dimCone->q[cntSubASoc] = numParameter + 1;
- cntSubASoc = cntSubASoc + 1;
- idxOptVarPhaseStart = idxOptVarPhaseStart + numPhaseOptVar[iphase];
- }
- socpInfo->cntNzAEqua = cntNzAEqua;
- socpInfo->cntRowAEqua = cntRowAEqua;
- socpInfo->cntNzAIneq = cntNzAIneq;
- socpInfo->cntRowAIneq = cntRowAIneq;
- socpInfo->cntNzASoc = cntNzASoc;
- socpInfo->cntRowASoc = cntRowASoc;
- socpInfo->cntSubASoc = cntSubASoc;
- }
- void getSocpInfo(cmscp_socpInfo** socpInfo0, cmscp_trajInfo** refTraj, cmscp_sparseInfo* sparseInfo, cmscp_consTypeInfo* consTypeInfo, cmscp_setup* setup)
- {
- // 函数功能:获取二阶锥规划问题的参数
- // 输入:
- // refTraj:参考轨迹信息
- // sparseInfo:稀疏(非零元位置)信息
- // consTypeInfo:约束类型信息
- // setup:问题设置参数
- // 输出:
- // socpInfo:存储有连续函数稀疏信息的结构体指针
- // 谢磊:2022 / 6 / 27编写
- cmscp_socpInfo* socpInfo;
- // 1. 获取二阶锥规划问题的维度信息和申请空间
- socpInfoSparseDim(socpInfo0, consTypeInfo, sparseInfo, setup);
- socpInfo = *socpInfo0;
- // 2. 获取二阶锥规划问题的目标函数
- socpInfoSparseObjEvent(socpInfo0, refTraj, sparseInfo, consTypeInfo, setup);
- // 3. 获取二阶锥规划问题的端点函数
- socpInfoSparseEndp(socpInfo0, refTraj, sparseInfo, consTypeInfo, setup);
- // 3. 获取二阶锥规划问题的连续函数
- socpInfoSparseCont(socpInfo0, refTraj, sparseInfo, consTypeInfo, setup);
- // 4. 获取二阶锥规划问题的衔接函数
- socpInfoSparseLink(socpInfo0, refTraj, sparseInfo, consTypeInfo, setup);
- // 5. 获取二阶锥规划问题的边界函数
- socpInfoBound(socpInfo0, consTypeInfo, setup);
- // 6. 获取二阶锥规划问题的信赖域函数
- socpInfoTrust(socpInfo0, refTraj, setup);
- // 7. 数据格式转化
- floatRcsToCcs(socpInfo->AEqua->ptrRow, socpInfo->AEqua->idxCol, socpInfo->AEqua->v,
- socpInfo->AEqua->m, socpInfo->AEqua->n, socpInfo->AEqua1->idxRow, socpInfo->AEqua1->ptrCol, socpInfo->AEqua1->v);
- floatRcsToCcs(socpInfo->ACone->ptrRow, socpInfo->ACone->idxCol, socpInfo->ACone->v,
- socpInfo->ACone->m, socpInfo->ACone->n, socpInfo->ACone1->idxRow, socpInfo->ACone1->ptrCol, socpInfo->ACone1->v);
-
- /* 打印二阶锥信息 */
- printIntVector(socpInfo->ACone->ptrRow, socpInfo->ACone->m + 1, "result\\AConePtrRow.txt");
- printIntVector(socpInfo->ACone->idxCol, socpInfo->ACone->ptrRow[socpInfo->ACone->m], "result\\AConeIdxCol.txt");
- printFloatVector(socpInfo->ACone->v, socpInfo->ACone->ptrRow[socpInfo->ACone->m], "result\\AConeValue.txt");
- printIntVector(socpInfo->AEqua->ptrRow, socpInfo->AEqua->m + 1, "result\\AEquaPtrRow.txt");
- printIntVector(socpInfo->AEqua->idxCol, socpInfo->AEqua->ptrRow[socpInfo->AEqua->m], "result\\AEquaIdxCol.txt");
- printFloatVector(socpInfo->AEqua->v, socpInfo->AEqua->ptrRow[socpInfo->AEqua->m], "result\\AEquaValue.txt");
- printFloatVector(socpInfo->bEqua, socpInfo->AEqua->m, "result\\bEqua.txt");
- printFloatVector(socpInfo->bCone, socpInfo->ACone->m, "result\\bCone.txt");
- printIntVector(socpInfo->dimCone->q, socpInfo->dimCone->nsoc, "result\\dimQ.txt");
- printFloatVector(socpInfo->c, socpInfo->numOptVar, "result\\c.txt");
- // 8 释放行压缩格式部分数据的空间
- freeSocpInfo(&socpInfo, setup, 0);
- }
|