LA 4255 Guess

本文探讨了一种解决特定拓扑排序问题的方法,通过构建有向图并进行拓扑排序,实现对给定数字序列进行正负判断。算法采用深度优先搜索策略,通过迭代调整节点值来寻找满足条件的解。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题目链接:https://icpcarchive.ecs.baylor.edu/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=2256

题意:给定一串数字a1,a2....an,给出sij的正负数值。sij代表ai+...aj的大小。求一组结果满足此条件。-10<=ai<=10.

思路:

设:Bi = a1+ a2 + ... ai,则,我们可以求出Bj和Bi-1的大小关系。可以由大到小连一条有向边,对此有向图做拓扑排序。

拓扑排序的做法如下:我们规定节点初始值:Bk =0(这只是一个相对值,事实上我们求的ai是Bi与Bi-1的差值,所以Bi初始是什么都没有关系,我们求得是相对差).

如果Bi 向Bj有一条有向边,那么Bi>Bj,说明Bi至少比Bj大1,我们用v[k]记录Bk的值,不断调整此数值。最终就能求得一组解。

#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <queue>
#include <iostream>

using namespace std;

#define Maxn 20
vector<int> g[Maxn];
int indeg[Maxn],v[Maxn];
int n;

void init()
{
	for(int i=0;i<Maxn;i++) 
	{
		g[i].clear();
		indeg[i] = 0;
		v[i] = 0;
	}
}

void toposort()
{
	queue<int> q;
	for(int i=0;i<=n;i++)
	{
		if(indeg[i] == 0) q.push(i);
	}
	while(!q.empty())
	{
		int s = q.front();
		q.pop();
		for(int i=0;i<g[s].size();i++)
		{
			int t = g[s][i];
			v[t] = min(v[t],v[s]-1);
			if(--indeg[t] == 0)
			{
				q.push(t);
			}
		}
	}
}
void output()
{
	for(int i=1;i<=n;i++)
	{
		if(i!=n) printf("%d ",v[i]-v[i-1]);
		else printf("%d\n",v[i]-v[i-1]);
	}
}
int main()
{
	#ifndef ONLINE_JUDGE
		freopen("in.txt","r",stdin);
	#endif
	int t;
	char str[20];
	scanf(" %d",&t);
	while(t--)
	{
		init();
		scanf(" %d",&n);
		scanf("%s",str);
		for(int i=1,p=0;i<=n;i++)
		{
			for(int j=i;j<=n;j++,p++)
			{
				if(str[p] == '-')
				{
					g[i-1].push_back(j);
					indeg[j]++; 
				}
				else if(str[p] == '+')
				{
					g[j].push_back(i-1);
					indeg[i-1]++;
				}
			}
		}
		toposort();
		output();
	}
	return 0;
}



clc; clear; close all; global setup % 基础参数设置(与原始代码相同) setup.n = [30]; % LG点个数 setup.ns = [5]; % 状态变量个数 setup.nu = [1]; % 控制变量个数 Re = 6378145; % 地球半径 rho0 = 1.225; % 海平面大气密度 H = 7200; % 大气范围 g0 = 9.80665; % 重力加速度 s = 15; % 参考面积 mach = 340; setup.Re = Re; setup.rho0 = rho0; setup.s = s; setup.H = H; setup.g0 = g0; setup.mach = mach; mI = 656080; % 运载器起飞质量 lI = 0.9; % 初始推重比 F0 = mI * lI * g0; F0_ = 1.1 * F0; setup.F0_ = F0_; % 过程约束 qmax = 50e3; nmax = 5; setup.qmax = qmax; setup.nmax = nmax; m0 = 136080; % 初始质量 mf = 60400; % 末端质量 t0 = 0; % 初始时间 t1 = 200; r0 = 2000 + Re; % 初始地心距 v0 = 200; % 初始速度 gamma0 = 30/57.3; % 初始弹道倾角 rf = 65000 + Re; % 终端高度 vf = 5 * mach; % 终端速度 gammaf = 20/57.3; % 终端弹道倾角 alpha0 = 20/57.3; % 起始攻角 alphaf = 0/57.3; % 终端攻角 % 无量纲化参数 tc = (Re/g0)^0.5; % 时间尺度 rc = Re; % 长度尺度 vc = (Re*g0)^0.5; % 速度尺度 Pc = m0 * g0; % 力尺度 mc = m0; % 质量尺度 setup.tc = tc; setup.rc = rc; setup.vc = vc; setup.Pc = Pc; setup.mc = mc; % 保存原始设置用于蒙特卡洛 orig_setup = setup; orig_m0 = m0; orig_r0 = r0; orig_v0 = v0; % 蒙特卡洛参数 num_runs = 100; % 运行次数 compute_times = zeros(num_runs, 1); % 存储计算时间 final_velocities = zeros(num_runs, 1); % 存储终端速度 % 优化选项设置 options = optimset('TolX',1e-20, 'TolFun',1e-20, 'TolCon',1e-20, ... 'MaxFunEvals',10000000, 'Algorithm','sqp', ... 'MaxIter',100, 'Display','off', 'largescale','on'); % 主循环 - 蒙特卡洛模拟 for run = 1:num_runs fprintf('蒙特卡洛运行 %d/%d...\n', run, num_runs); % 重置全局设置 setup = orig_setup; % 添加随机扰动 (±5%) m0_perturbed = orig_m0 * (0.95 + 0.1*rand()); r0_perturbed = orig_r0 * (0.95 + 0.1*rand()); v0_perturbed = orig_v0 * (0.95 + 0.1*rand()); % 重新计算无量纲初始状态 r0dot = r0_perturbed / rc; v0dot = v0_perturbed / vc; gamma0dot = gamma0; m0dot = m0_perturbed / mc; alpha0dot = alpha0; X0 = [r0dot; v0dot; gamma0dot; m0dot; alpha0dot]; setup.X0 = X0; setup.r0dot = r0dot; setup.v0dot = v0dot; setup.m0dot = m0dot; setup.gamma0dot = gamma0dot; setup.alpha0dot = alpha0dot; % 重新初始化阶段参数 phase = struct(); iphase = 1; phase.t0{iphase} = t0/tc; phase.tf{iphase} = t1/tc; phase.m0{iphase} = m0_perturbed/mc; phase.mf{iphase} = mf/mc; % 生成Legendre-Gauss点 [phase.tau{iphase}, phase.w{iphase}, phase.D{iphase}] = LG_Points(setup.n(iphase)); % 初始猜测 phase.X{iphase} = [r0dot v0dot gamma0dot phase.m0{iphase} alpha0dot; r0dot v0dot gamma0dot phase.mf{iphase} alpha0dot]; phase.Xguess{iphase} = interp1([-1;1], phase.X{iphase}, [-1; phase.tau{iphase}; 1], 'spline'); phase.Uguess{iphase} = (1/57.3) * ones(setup.n(iphase), 1); phase.x0guess{iphase} = [reshape(phase.Xguess{iphase}, [], 1); reshape(phase.Uguess{iphase}, [], 1)]; % 边界约束 rmindot = Re / rc; rmaxdot = 2 * Re / rc; vmindot = -10000 / vc; vmaxdot = -vmindot; gammamindot = 0; gammamaxdot = 80/57.3; alphamindot = -10/57.3; alphamaxdot = 30/57.3; dalphamin = -5/57.3; dalphamax = 5/57.3; phase.XU{iphase} = [rmaxdot*ones(setup.n(iphase)+2,1) vmaxdot*ones(setup.n(iphase)+2,1) ... gammamaxdot*ones(setup.n(iphase)+2,1) phase.m0{iphase}*ones(setup.n(iphase)+2,1) ... alphamaxdot*ones(setup.n(iphase)+2,1)]; phase.UU{iphase} = dalphamax * ones(setup.n(iphase), 1); phase.xU{iphase} = [reshape(phase.XU{iphase}, [], 1); reshape(phase.UU{iphase}, [], 1)]; phase.XL{iphase} = [rmindot*ones(setup.n(iphase)+2,1) vmindot*ones(setup.n(iphase)+2,1) ... gammamindot*ones(setup.n(iphase)+2,1) phase.mf{iphase}*ones(setup.n(iphase)+2,1) ... alphamindot*ones(setup.n(iphase)+2,1)]; phase.UL{iphase} = dalphamin * ones(setup.n(iphase), 1); phase.xL{iphase} = [reshape(phase.XL{iphase}, [], 1); reshape(phase.UL{iphase}, [], 1)]; setup.phase = phase; x0 = [phase.x0guess{1}; phase.tf{1}]; xU = [phase.xU{1}; phase.tf{1}]; xL = [phase.xL{1}; phase.t0{1}]; % 运行优化并计时 tic; [x, f] = fmincon(@Objective, x0, [], [], [], [], xL, xU, @NonlconFun, options); compute_times(run) = toc; % 提取终端速度(无量纲转有量纲) final_velocities(run) = x(2) * vc; fprintf('运行 %d 完成: 时间 = %.2f 秒, 终端速度 = %.2f m/s\n', ... run, compute_times(run), final_velocities(run)); end % 生成散点图 figure('Color', 'white', 'Position', [100, 100, 800, 600]); scatter(compute_times, final_velocities, 50, 'filled', 'MarkerFaceColor', [0.2, 0.4, 0.8], ... 'MarkerEdgeColor', 'k', 'LineWidth', 1.2); % 添加趋势线 hold on; p = polyfit(compute_times, final_velocities, 1); yfit = polyval(p, compute_times); plot(compute_times, yfit, 'r-', 'LineWidth', 2); % 图表标注 title('轨迹优化速度性能分析', 'FontSize', 16, 'FontWeight', 'bold'); xlabel('计算时间 (秒)', 'FontSize', 14); ylabel('终端速度 (m/s)', 'FontSize', 14); grid on; box on; % 添加统计信息 mean_time = mean(compute_times); mean_vel = mean(final_velocities); text(0.7*max(compute_times), min(final_velocities)+100, ... sprintf('平均时间: %.2f s\n平均速度: %.2f m/s', mean_time, mean_vel), ... 'FontSize', 12, 'BackgroundColor', 'white'); % 保存结果 save('monte_carlo_results.mat', 'compute_times', 'final_velocities'); disp('蒙特卡洛模拟完成,结果已保存。'); % 辅助函数:生成Legendre-Gauss点 function [tau, w, D] = LG_Points(n) [tau, w] = legendre_points(n); D = collocation_matrix(tau); end function [x, w] = legendre_points(n) % 返回n阶Legendre多项式的根(LG点)和权重 beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2)); % 三项递推系数 T = diag(beta,1) + diag(beta,-1); % Jacobi矩阵 [V, D] = eig(T); % 特征值分解 x = diag(D); % LG点 w = 2*V(1,:).^2; % 权重 % 排序 [x, idx] = sort(x); w = w(idx)'; end function D = collocation_matrix(tau) % 计算微分矩阵 n = length(tau); D = zeros(n); for j = 1:n for k = 1:n if j ~= k D(j,k) = lagrange_derivative(tau, j, k); end end end D(j,j) = -sum(D(j,:)); end function d = lagrange_derivative(tau, j, k) % 计算Lagrange基函数在节点处的导数 n = length(tau); prod_val = 1; for m = 1:n if m ~= j && m ~= k prod_val = prod_val * (tau(k) - tau(m))/(tau(j) - tau(m)); end end d = prod_val / (tau(j) - tau(k)); end % 目标函数 (简化示例) function J=Objective(x) global setup J=x(end)*setup.tc; end % 非线性约束函数 (简化示例) function [c, ceq] = NonlconFun(x) % 实际应用中应根据具体问题实现 c = []; % 无不等式约束 ceq = []; % 无等式约束 end在以上代码中加入以下非线性约束函数:%% 约束满足 function [c,ceq]=NonlconFun(x) global setup rho0=setup.rho0; Re=setup.Re; H=setup.H; s=setup.s; g0=setup.g0; qmax=setup.qmax; nmax=setup.nmax; %无量纲化参数 rc=setup.rc; vc=setup.vc; tc=setup.tc; mc=setup.mc; n=setup.n; %LG点个数 ns=setup.ns; %状态变量个数 nu=setup.nu; %控制变量个数 np=length(n); %阶段数 nx=[0]; nx(1)= n(1)*(ns(1)+nu(1))+2*ns(1); setup.nx=nx; xi{1}=x(1:nx(1)); Xi{np}=[]; Ui{np}=[]; Fi{np}=[]; DXi{np}=[]; t0=setup.phase.t0; tf=setup.phase.tf; D=setup.phase.D; w=setup.phase.w; oe=setup.oe; tf{np}=x(end); if tf{np}<1/tc aa=0; end for i1=1:np X=ones(n(i1)+2,ns(i1)); for i2=1:ns(i1) X(:,i2)=xi{i1}((i2-1)*(n(i1)+2)+1:i2*(n(i1)+2),1); end Xi{i1}=X; U=ones(n(i1),nu(i1)); for i3=1:nu(i1) U(:,i3)=xi{i1}(ns(i1)*(n(i1)+2)+(i3-1)*n(i1)+1:ns(i1)*(n(i1)+2)+i3*n(i1),1); end Ui{i1}=U; F=MyDeo(X,U,i1); Fi{i1}=F; %动力学约束 ceqd=D{i1}*X(1:end-1,:)-(tf{i1}-t0{i1})/2*F; ceqdr=reshape(ceqd,[],1); ceqdri{i1}=ceqdr; %终端约束 ceqe=X(1,:)+(tf{i1}-t0{i1})/2*w{i1}'*F-X(end,:); ceqei{i1}=ceqe'; %路径约束 v=X(2:end-1,2); %动压 r=X(2:end-1,1); rho=rho0*exp(-(r-1)*Re./H); q=0.5.*rho.*(v.*vc).^2; cq=q-qmax*ones(n(i1),1); ciq=cq; ceqi{i1}=[ceqdr;ceqe']; ci{i1}=[-1]; end ceql=[0]; %边界条件-初值 ceqs=Xi{1}(1,:)-setup.X0'; ceqsr=ceqs'; %边界条件-终值 rf=Xi{1}(end,1)'; Vf=Xi{1}(end,2)'; ceqo=[Vf]-oe(2); ceq=[ceqi{1};ceqsr;ceql;ceqo']; c=[ci{1}]; end 生成完整的代码
最新发布
06-13
<think>我们正在处理一个优化问题,需要在MATLAB中集成非线性约束函数NonlconFun。根据引用[1]和[2],在MATLAB中求解非线性优化问题时,可以通过fmincon函数来添加非线性约束。非线性约束函数需要返回两个输出:不等式约束c和等式约束ceq。同时,我们需要确保非线性约束函数的输入为优化变量x,并且返回的c和ceq满足c<=0和ceq=0。假设我们已有的代码框架如下(使用fmincon求解优化问题):x0=...;%初始值A=...;%线性不等式约束矩阵b=...;%线性不等式约束向量Aeq=...;%线性等式约束矩阵beq=...;%线性等式约束向量lb=...;%变量下界ub=...;%变量上界options=optimoptions('fmincon','Display','iter','Algorithm','sqp');%设置选项[x,fval]=fmincon(@(x)ObjectiveFun(x),x0,A,b,Aeq,beq,lb,ub,@(x)NonlconFun(x),options);其中,非线性约束函数NonlconFun需要单独编写,其基本格式为:function[c,ceq]=NonlconFun(x)%计算非线性不等式约束,c<=0c=...;%可以是多个约束组成的向量%计算非线性等式约束,ceq=0ceq=...;%可以是多个约束组成的向量end现在,我们需要将用户给定的非线性约束函数NonlconFun集成到上述框架中。步骤:1.确保NonlconFun.m文件与主脚本在同一目录下,或者将其放在MATLAB的搜索路径中。2.在主优化代码中,通过函数句柄@NonlconFun来引用这个函数。3.确保NonlconFun函数接受一个向量输入x,并返回两个向量:c(不等式约束)和ceq(等式约束)。注意:如果约束函数需要额外的参数,可以使用匿名函数来传递,例如:[x,fval]=fmincon(@(x)ObjectiveFun(x,param1,param2),...,@(x)NonlconFun(x,param1,param2),...);下面我们给出一个完整的示例代码,包括目标函数和非线性约束函数的定义,以及调用fmincon的过程。示例:假设我们要求解以下优化问题:目标函数:$f(x)=\exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1)$约束条件:非线性不等式约束:$1.5+x(1)*x(2)-x(1)-x(2)\leq0$$-x(1)*x(2)\leq10$非线性等式约束:$x(1)^2+x(1)*x(2)-10=0$则非线性约束函数NonlconFun.m可以写为:function[c,ceq]=NonlconFun(x)%不等式约束(注意要写成<=0的形式)c=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];%第二个约束改写为:-x1*x2-10<=0等价于-x1*x2<=10%等式约束ceq=x(1)^2+x(1)*x(2)-10;end目标函数文件ObjectiveFun.m:functionf=ObjectiveFun(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);end主脚本(例如main.m):x0=[-1,1];%初始值A=[];%无线性不等式约束b=[];Aeq=[];%无线性等式约束beq=[];lb=[];%无下界ub=[];%无上界%设置优化选项options=optimoptions('fmincon','Display','iter','Algorithm','interior-point');%调用fmincon[x,fval]=fmincon(@ObjectiveFun,x0,A,b,Aeq,beq,lb,ub,@NonlconFun,options);运行主脚本即可得到优化结果。如果用户需要将非线性约束函数直接写在主脚本中(作为局部函数),则可以将NonlconFun定义在主脚本的末尾(在MATLABR2016b及以上版本支持局部函数)。例如,在同一个m文件中:%主脚本部分x0=[-1,1];...%同上[x,fval]=fmincon(@objfun,x0,A,b,Aeq,beq,lb,ub,@nonlcon,options);%目标函数(局部函数)functionf=objfun(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);end%非线性约束函数(局部函数)function[c,ceq]=nonlcon(x)c=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];ceq=x(1)^2+x(1)*x(2)-10;end注意事项:1.确保约束函数的输出c和ceq为列向量(如果是多个约束)。2.如果某个约束不存在,则返回空数组,例如:ceq=[]。3.非线性约束函数中不要忘记将约束写成标准形式(c<=0,ceq=0)。4.使用合适的算法(如'interior-point','sqp'等),有些算法对非线性约束有更好的支持。如果用户已有目标函数和约束函数的表达式,按照上述步骤即可集成。针对用户的问题,如果已经有一个优化框架,只需要将NonlconFun函数按照要求编写,并在fmincon调用时作为第9个参数传入即可。如果用户需要完整的代码,请提供目标函数、非线性约束的具体表达式以及线性约束、边界条件等,我们可以帮助编写。</think>### 在MATLAB中集成非线性约束函数的完整指南 要将非线性约束函数 `NonlconFun` 正确集成到MATLAB优化代码中,请遵循以下步骤: #### 1. 非线性约束函数格式要求 非线性约束函数必须满足MATLAB标准格式: ```matlab function [c, ceq] = NonlconFun(x) % 非线性不等式约束 (c(x) ≤ 0) c = []; % 非线性等式约束 (ceq(x) = 0) ceq = []; end ``` - `c`:不等式约束向量,需满足 $c(x) \leq 0$ - `ceq`:等式约束向量,需满足 $ceq(x) = 0$ - `x`:优化变量向量 #### 2. 完整优化代码框架 ```matlab %% 优化问题设置 x0 = [1, 2]; % 初始点 A = []; % 线性不等式系数矩阵 b = []; % 线性不等式右侧向量 Aeq = []; % 线性等式系数矩阵 beq = []; % 线性等式右侧向量 lb = [0, 0]; % 变量下界 ub = [10, 10]; % 变量上界 %% 目标函数定义 objective = @(x) x(1)^2 + x(2)^2; % 示例目标函数 $f(x) = x_1^2 + x_2^2$ %% 优化选项配置 options = optimoptions('fmincon', ... 'Display', 'iter', ... % 显示迭代过程 'Algorithm', 'sqp', ... % 选择序列二次规划算法 'MaxIterations', 1000, ... % 最大迭代次数 'StepTolerance', 1e-6); % 步长容差 %% 执行优化(集成非线性约束) [x_opt, fval, exitflag, output] = fmincon(... objective, ... % 目标函数 x0, ... % 初始点 A, b, ... % 线性不等式约束 Aeq, beq, ... % 线性等式约束 lb, ub, ... % 变量边界 @NonlconFun, ... % 非线性约束函数 options); % 优化选项 %% 结果输出 disp('最优解:'); disp(x_opt); disp(['最优目标值: ', num2str(fval)]); disp(['退出状态: ', num2str(exitflag)]); ``` #### 3. 非线性约束函数示例 假设需要满足约束: 1. $x_1x_2 - 20 \leq 0$ (非线性不等式) 2. $x_1^2 + x_2 - 5 = 0$ (非线性等式) 对应的 `NonlconFun.m`: ```matlab function [c, ceq] = NonlconFun(x) % 非线性不等式约束 (c ≤ 0) c = x(1)*x(2) - 20; % $c(x) = x_1x_2 - 20 \leq 0$ % 非线性等式约束 (ceq = 0) ceq = x(1)^2 + x(2) - 5; % $ceq(x) = x_1^2 + x_2 - 5 = 0$ end ``` #### 4. 关键注意事项 1. **约束标准化**:确保所有约束转换为标准形式: - 不等式:$g(x) \leq 0$ - 等式:$h(x) = 0$ 2. **梯度计算**(可选): 对于复杂问题,提供解析梯度可加速收敛: ```matlab function [c, ceq, gradc, gradceq] = NonlconFun(x) c = x(1)*x(2) - 20; ceq = x(1)^2 + x(2) - 5; % 不等式约束梯度 gradc = [x(2); x(1)]; % $\nabla c(x) = [x_2, x_1]^T$ % 等式约束梯度 gradceq = [2*x(1); 1]; % $\nabla ceq(x) = [2x_1, 1]^T$ end ``` 在选项中启用梯度: ```matlab options = optimoptions('fmincon', 'SpecifyConstraintGradient', true); ``` 3. **错误排查**: - 检查约束函数返回值是否为实数 - 确保约束函数在定义域内处处可计算 - 使用 `try-catch` 块捕获异常: ```matlab try [x_opt, fval] = fmincon(...); catch ME disp('优化错误:'); disp(ME.message); % 添加调试断点 keyboard; end ``` #### 5. 完整工作流验证 1. 创建目标函数文件 `objective.m` 2. 创建非线性约束文件 `NonlconFun.m` 3. 主脚本集成优化框架 4. 运行脚本并检查输出: - `exitflag > 0` 表示成功收敛 - 验证解是否满足所有约束 > 非线性约束的集成是优化问题求解的关键环节,合理定义约束函数可显著提高求解效率和精度[^4]。对于复杂约束系统,建议采用分段调试策略,先验证单个约束再组合测试[^3]。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值